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RESEARCH SUMMARY 


A simulator for predicting even- and uneven-aged stand 
development was constructed and validated for the ponderosa 
pine/Arizona fescue habitat type of the Southwest. The structure 
and dynamics of the simulator make it potentially useful for 
answering both even- and uneven-aged management questions. 


The primary framework for the simulator is two diameter 
distributions, one for blackjack pine and the other for yellow 
pine (both are different vigor classes of ponderosa pine). The 
diameter distributions are expressed as number of trees in each 
l-inch diameter class, and the lowest diameter class is 4 inches 
for uneven-aged stands and 1 inch for even-aged stands. Within 
this framework, there are four basic components: upgrowth, 
mortality, conversion, and ingrowth. i 


Upgrowth is the movement of trees from a given diameter 
class into larger diameter classes. It was modeled as a function 
of the distribution of trees within the diameter class and of 
basal area growth of the diameter class. During modeling, diam- 
eter class basal area growth rates of even- and uneven-aged 
stands were shown to be logically interrelated if the structure 
of within-stand competition is incorporated in the model. 


Diameter class mortality was modeled as the proportion of 
trees dying to eliminate the problem of negative survival rates. 
Even- and uneven-aged endemic mortality rates also proved to be 
logically interrelated if catastrophic mortality was treated in a 
fashion similar to cutting "mortality." 


Conversion is the process in which blackjack pines transform 
into yellow pines. As with mortality, conversion was modeled as 
a proportion of the diameter class transforming. 


The final component of the simulator is ingrowth, which is 
the number of trees growing into the 4- or 5-inch diameter classes 
of blackjack pine. Data restrictions precluded the development 
of a regeneration model and limited the ingrowth model to uneven- 
aged stands. For even-aged stands, it is assumed that no (or 
insignificant) trees exist below the l-inch diameter class. 
Simulator applications are further restricted to stands with an 
average stand diameter of 4 inches or greater. 


Two equations were needed to model the ingrowth process. 
The first predicts the total number of trees ingrowing into the 
4- and 5-inch diameter classes. The second equation predicts 
the proportion of total ingrowth that will through grow into the 
5-inch diameter class in the same growth period. 


During the validation process, an assumption made by past 
uneven-aged modelers was tested: that small plots have the same 
structure and dynamics as the total stand and can therefore be 
used to develop whole stand simulators. From this analysis, I 
concluded that developing a simulator, using 2.5-acre plots, and 
then using it to predict stand averages generally reduces both 
the precision and accuracy of the stand estimates. In this 
study, however, the loss appears not too severe for most potential 
uses of a whole stand simulator. 
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STUDY OBJECTIVES 


Recent publications (USDA Forest Service 1975, 1976; Alexander and Edminster 1977a, 1977b; 
Hann and Bare 1979) suggest a renewal of interest in uneven-aged forest management. Hann and 
Bare (1979) reviewed major questions facing a forest manager interested in practicing uneven- 
aged management and summarized the analytical tools currently available to answer these ques- 
tions. They demonstrated the important role stand development simulators play in making 
decisions using some of the more recent analytical tools. The construction of appropriate 
stand development simulators is therefore a prerequisite to the application of these analytical 
tools. 


A number of simulator types have been developed in recent years (Munro 1974). But techno- 
logical limitations of the available uneven-aged analytical tools require using whole-stand 
simulators that characterize the stand by at least the number of trees in various diameter 
classes (Adams 1974; Adams and Ek 1974, 1976; Hann 1978). With this constraint, the objectives 
of this study became: 


1. To examine the dynamics of uneven-aged stand development in the ponderosa pine/ 
Arizona fescue habitat type by developing a whole-stand simulator that would have the potential 
for aiding forest managers in answering questions concerning uneven-aged management. Ponderosa 
pine was chosen because it is one of.the few species in the West for which permanent plot data 
exist for both even- and uneven-aged stands. 


2. To examine whether even- and uneven-aged stand development in ponderosa pine/Arizona 
fescue could logically be interrelated and, if so, to develop the whole-stand simulator to be 
applicable to both. 


3. To validate the whole-stand simulator by analyzing its capability to predict both 
average plot and average stand development. 


PREVIOUS UNEVEN-AGED, WHOLE-STAND SIMULATORS 


Two approaches to the development of uneven-aged, whole-stand simulators have been 
reported. The first approach uses differential equations and the second uses difference 
equations to characterize stand dynamics. The differential equation approach is based on using 
instantaneous rate equations to model stand dynamics. Two alternate techniques exist for 
parameter estimation of the instantaneous rate equations. In the first technique, modelers 
approximate the instantaneous rates through use of periodic rates that have been converted 
to an annual basis, and then estimate the parameters through regression analysis (Moser and 
Hall 1969; Moser 1972, 1974). If the length of the growth period is short and if the instan- 
taneous rates do not change drastically within the period, then this technique can provide 
satisfactory estimates of the instantaneous rate equations. The second technique uses the 
multipoint boundary value method to develop parameter estimates of the instantaneous rate 
equations (Leary 1970). This process is a common numerical analysis method but requires data 
from at least three different points in the development of the stand. While Leary's (1970) use 
of this method was in even-aged stands, it would also be appropriate in uneven-aged stands. 


The second approach recognizes that most forest rate data are periodic in nature and, 
therefore, an appropriate way to model them is by using difference equations. An example of 
an uneven-aged, whole-stand simulator using this approach is Ek's (1974) mixed species, north- 
ern hardwood simulator that characterizes the stand by number of trees in 2-inch diameter 
classes. In this simulator, the number of trees in each diameter class at the end of a growth 
period is a linear combination of ingrowth, upgrowth, and mortality. Ingrowth is the number 
of trees growing into the smallest diameter class in a growth period, and Ek modeled it as a 
nonlinear function of total stand basal area and total number of trees. Upgrowth is the 
number of trees growing out of a specified diameter class in a growth period, and it was 
modeled as a nonlinear function of number of trees and basal area in the diameter class, total 
number of trees and basal area in the stand, and site index. Mortality is the number of trees 
in a specified diameter class that die in a growth period, and it was modeled as a nonlinear 
function of number of trees and basal area in the diameter class, and total number of trees 
and basal area in the stand. Equations were also developed for predicting rough cordwood 
volume and board-foot volume. Cuttings were introduced by reducing the number of trees in the 
appropriate diameter class at the end of the growth period. 


As a demonstration of the application of the uneven-aged analytical decisionmaking tools 
they had developed, Adams and Ek (1974) used a modified version of Ek's (1974) model. 


The difference equation approach was used in this study for two reasons. First, the 
development of difference equations appeared to be the most logical use of the periodic data 
available. Because the data source and the equation type are both periodic, it is not necessary 
to make numerical approximations in order to develop and/or apply the equations. Most models 
using differential equations do require numerical approximation methods to estimate parameters 
and/or to solve them. 


Second, the use of difference equations avoids potential problems of inappropriately 
applying the equations to growth periods different from that used in equation development. 
Because of the continuity of differential equations, it is possible to predict growth for any 
period. The differential equations are usually, however, approximations of the true continuous 
growth process (not only can the parameter estimates be approximations, but the model forms 
themselves are often approximations). As such, they should not be applied to growth periods 
much different from those used in equation development or else erroneous results can occur. On 
the other hand, in normal use difference equations can only be used to predict growth for 
periods equal to those used to develop the equations. Multiple periodic growth (or fractional 
periodic growth, if the user wishes to chance extending the equations in that fashion) is 
derived from these periodic growth rates. 


STUDY DATA BASE 


Plot Characteristics 


The data used in this study were from the Fort Valley Experimental Forest, located a few 
miles northwest of Flagstaff, Arizona. The uneven-aged data consist of long-term remeasure- 
ments on four 80-acre plots. Each plot has been divided into 32 or 34 subplots ranging in size 
from 0.5 to 3.1 acres, but most of these Subplots average near 205 acres in size. ~ The measune— 
ment of the plots started either in 1920 or 1925, and remeasurements have occurred at 5-year 
intervals until the 1960's when the interval was lengthened to 10 years. The two attributes 
recorded on every tree above the lower diameter limit were diameter at breast height (d.b.h.) 
and tree condition (a classification as to presence of a damaging agent if the tree is alive, 
or to the cause of mortality if the tree has died within the previous growth period). Age 
class and age-vigor class of each tree were measured every fourth growth period. Age class 
categorizes a tree as being blackjack or yellow pine, and age-vigor uses age and size infor- 
mation, in combination with the previous 20-year average diameter growth rate, to help identify 
the tree's growth potential. The lower diameter limit was first set at 3.6 inches (the lower 
bound for the 4-inch diameter class). It was then changed to 7.6 inches in 1940, and recently 
to 6.0 inches. A brief history, and description’ of “cach uneven-agedy plot 1s im stable n: 


Table 1.--History and description of uneven-aged study plots 


Plot 

No. Description 

61 Represents virgin ponderosa pine. The south half of the plot was first 
measured in 1920 and the other half in 1925. Average site index of this plot 
has been estimated to be approximately 67 by Minor's (1964) site curve. 

62 First measured in 1920. It was salvage cut in 1940 and again in 1953. 
Average site index is 73 (Minor 1964). 

vat Had a modified selection cut in 1924, followed by its first measurement in 1925. 
A timber stand improvement (T.S.I.) cut occurred in 1935, an improvement 
selection cut in 1946, a precommercial thinning in 1967, and a group selection 
cut in 1968. The basal area stocking level now established as a goal for this 
plot is 96 square feet. Average site index is 73 (Minor 1964). 

TZ. Treatment has been identical to that of plot 71. The basal area stocking level 


now set for this plot, however, is 70 square feet. Average site index is 70 
(Minor 1964). 


The even-aged data came from the Taylor Woods level of growing stock study, located very 
close to the uneven-aged plots (Schubert 1971), and consisting of 18 plots ranging in size from 
0.75 to 1.25 acres. These plots represent three replications of six growing stock levels (30, 
60, 80, 100, 120, and 150 square feet of basal area). The plots were measured in 1962, 1967, 
and 1972. They were thinned in 1962 and 1972 (the 1962 thinning was before measurement). All 
trees on the plots were measured. Average site index for the study area is 88 (Minor 1964). 


The attributes measured on each tree include d.b.h., tree condition, crown class (whether 
a tree was a dominant, codominant, and so forth), crown quality (the condition of the crown), 
and percent of the bole in live crown. On a subsample of the trees, the following additional 
items were measured: total tree height, crown length, and crown width. 


Data Problems 


While this data source is unusual for the West because it provides both even-aged data and 
long-term uneven-aged growth data, it also posed numerous modeling problems. A brief discus- 
sion follows on these problems and the actions taken to minimize them. See appendix A for a 
more complete discussion. 


1. All of the data used in this study came from a narrow geographic and site range. 
Therefore, it is recommended that the simulator developed in this study be applied only to 
stands on the ponderosa pine/Arizona fescue (Pinus ponderosa/Festuca arizonica) habitat type of 
northern Arizona. 


2. Because of changes in the lower limit of diameter measurements, only data measured up 
to 1940 were used to estimate component model parameters. 


3. An intrusion of a major paved highway through all of the uneven-aged plots necessi- 
tated the removal of some subplots. 


4. As mentioned earlier, a wide range of subplot sizes exists on the uneven-aged plots. 
To avoid a potential problem of confounding subplot stand differences with subplot size differ- 
ences, all uneven-aged subplots below 2 acres were eliminated from the data base. 


5. Because of a lack of height-growth data on the uneven-aged plots, development of a 
height-growth model was not possible. Instead, additional data needed for the development of 
height-diameter equations were obtained from the School of Forestry at Northern Arizona Uni- 
versity and from the Rocky Mountain Forest and Range Experiment Station at Flagstaff. 


Plot and Subplot Allocation 


With the elimination of the small subplots and those seriously affected by the highway, 
110 uneven-aged subplots remained. Of these, 24 were reserved for validation (objective 
number 3). To pick these, each of the four uneven-aged plots had their subplots divided into 
low, medium, and high basal area stocking classes (one-third of the subplots in each), and two 
subplots were randomly selected from each class. This helps to insure that the validation data 
represent the full range of stocking conditions. 


Also reserved for validation were six even-aged plots. These were also randomly chosen by 
selecting one plot from each of the six basal area stocking classes. The result was 86 uneven- 
aged and 12 even-aged plots for use in developing the simulator. 


SIMULATOR DEVELOPMENT 


The development of a whole-stand simulator that can express number of trees by diameter 
class is the main objective of this study. A l-inch diameter class size was chosen because it 
offered a good compromise between model complexity on one hand and sensitivity to stand struc- 
ture and to managerial and silvicultural information needs on the other. The decision to 
express the stand as number of trees instead of basal area in each diameter class was made 
because many of the simulator's stand dynamic components are best expressed in terms of 
number of trees. 


The dynamics of stand development is composed of three basic elements: growth, mortality, 
and recruitment. 


Growth 


The process of growth takes place on all portions of the tree: roots, crown, and stem. 
The growth of most interest to the manager, however, is stem or bole growth, which is the 
height and diameter growth. With the lack of height-growth information in the data set, a 
height-growth model could not be developed. Instead, height equations were developed as a 
substitute to aid in the prediction of product potential. Lack of a height-growth model did 
not affect the processes used to predict change in number of trees by diameter class. Such a 
change is the difference between the number of trees growing into the diameter class and the 
number growing out of it. The number of trees growing from one diameter class to another, 
sometimes called upgrowth, is a function of three components: diameter growth, the distribu- 
tion of number of trees within the diameter class, and diameter class size (Wahlenberg 1941). 


At least two approaches exist to modeling this process. In one approach, Ek (1974) chose 
not to treat these components separately, but rather, he predicted the number of trees growing 
into the next larger diameter class directly. Another approach is to consider each component 
of the upgrowth model individually. This approach was used in the study because it allowed a 
closer examination of the validity of the individual component assumptions incorporated into 
the upgrowth model. 


DISTRIBUTION WITHIN DIAMETER CLASS AND DIAMETER CLASS SIZE 


The distribution of the trees within a diameter class can influence the number of trees 
advancing out of the class, given an average class diameter growth rate (Wahlenberg 1941; Meyer 
1942). The easiest distribution to use, and the one assumed by Moser (1974), is the uniform 
distribution. If this distribution is valid, then the following method can be used for comput- 
ing upgrowth (Wahlenberg 1941; Meyer 1942). The quotient found by dividing diameter growth by 
diameter class size is computed. The integer portion of this quotient is then the number of 
diameter classes in which all trees will advance, and the fractional portion of the quotient is 
the proportion of the trees that will advance one diameter class further. (For example, if the 
quotient was 2.29, then all trees will advance two diameter classes, and 29 percent of them 
will advance three diameter classes.) 


Because it is the easiest to apply, I first tested the appropriateness of the uniform 
distribution through use of the chi-square "goodness-of-fit" test (see appendix B for details). 
A total of 714 diameter classes were tested at the 99 percent level of significance, and of 
these, only seven were found to be significantly different from the uniform distribution. 

These findings were further substantiated by visually checking the distributions of many of 
these classes. Therefore, the assumption of within-class uniformity could not be rejected in 
the even- and uneven-aged stands of the Fort Valley Experimental Forest. 


DIAMETER GROWTH 


The final component necessary for the calculation of upgrowth is an estimate of the diam- 
eter growth rate for each class. The following summarizes procedures used to develop the 
appropriate diameter growth models. A more complete discussion is in appendix C. 


Initial Review and Modeling Decisions 


A review of the literature for southwestern ponderosa pine helped identify factors histor- 
ically recognized as related to diameter growth. The significant factors identified included: 
site index, rainfall, total stand density, diameter class size, tree vigor as indicated by bark 
color (in other words, "yellow" pine and "blackjack" pine), and the structure of the competition 
as indicated by the position of the diameter class in the stand. 


Given the general factors that can affect diameter growth, it was then necessary to 
choose: the general structure of the model, the form of the dependent variable, the method for 
handling the error structure of the model during prediction, and the appropriate transforma- 
tions of the independent variables. After careful consideration, I decided to predict the 
natural log of diameter class basal area growth and then convert that value to a diameter 
growth rate. Three basic reasons for this choice are: 


1. A nonlinear, multiplicative error model form probably best represented the interaction 
of the independent variables with themselves and their effect upon diameter growth. 


2. Basal area growth is often nearly linear over short time periods, which simplifies 
the extrapolation of growth rates to growth periods different from that originally used in 
equation development. 

3. The residuals of log of basal area growth more often approach normality and homo- 
geneous variance than other dependent variables, which is advantageous in model development and 
significance testing. 


Several ways of using information about the error structure of the model are available when 
using the finished model for predictive purposes. In deterministic modeling, the random error 
element of the model is simply ignored for predictive purposes. If the model is properly 
constructed, predictions from a deterministic model represent the "expected" value of the 
Output. In stochastic modeling, the random error element is considered in prediction. To 
obtain an expected value of the output from a stochastic model, however, requires a ''Monte 
Carlo analysis," which can be time consuming and expensive. 


Stage (1973) described an alternative method, which could be termed a ''Monte Carlo swindle," 
of introducing a stochastic element into his individual tree log of basal area growth equations 
while maintaining the relative simplicity of a deterministic model. Because it was expected 
that individual tree growth within a diameter class would be quite variable, elements from 
Stage's (1973) "Monte Carlo swindle'' were adapted and used in this study. 


The resulting process consists of dividing the number of trees in each diameter class into 
thirds (one each to represent the fast, slow, and "moderate" srowers)),~predicting the! expected 
residual for each third, and adding or subtracting this to the predicted average log of basal 
area growth value for the diameter class. In this fashion, estimates of basal area growth are 
obtained for each third of the diameter class. Upgrowth is then calculated separately for each 
third, which allows the trees in a diameter class to move into a wider range of larger diameter 
classes. More realistic predictions should result. 


The final decision before parameter estimation could begin was the choice of appropriate 
transformations necessary to create those independent variables most likely to be highly corre- 
lated to the dependent variable. Information used in selecting the independent variables 
included the silvicultural literature on diameter growth factors for southwestern ponderosa 
pine and prior general ''growth and yield' knowledge reported in the mensurational literature. 
From this work, 43 potentially useful independent variables were identified and examined in the 
next phase of the diameter growth modeling process. 


Development of Equations 


Ordinary, least squares regression techniques were used to develop separate equations for 
both blackjack and yellow pine. I initially separated these two vigor classes because evidence 
from prior studies suggest significant growth differences between blackjack and yellow pine. 
The analysis process to determine the form and parameters of the finished models required six 
phases: 


Phascale An all-combinations screening run was made On each individual plot 
data set to eliminate those independent variables not highly corre- 
lated to the dependent variable. 

Phase 2. The separate plot data were then combined into three sets: virgin, 
uneven-aged data; managed, uneven-aged data; and managed, even- and 
uneven-aged data. A second set of all combination screening runs 
was then made for each data set and for yellow and blackjack pine. 
From these runs, a common set of independent variables exhibiting 
behavior both reasonable and consistent with mensurational and 
silvicultural expectations was selected for both yellow and black- 
jack pine. Because the basic model forms and the size.of the 
parameter estimates differed so greatly between yellow and black- 
jack pine, I decided that at this phase of the analysis any 
attempt to combine the two vigor classes through use of dummy 
variables would be futile. Therefore, separate equations were 
maintained throughout the rest of the analysis. 

Phases. The combined data sets were further collapsed into two sets: uneven- 
aged and even- and uneven-aged. The reason for maintaining separate 
sets was to allow comparison of the effect upon predictive capability 
of adding the even-aged data to the uneven-aged. To do this merging, 
the differences between the model coefficients for the managed and 
the virgin data sets were handled as functions of time-since-last- 
cutting. Sigmoidal transformations of time-since-last-cutting 
were used because the effect of cutting was expected to asymptotically 
approach the virgin growth rate as time-since-last-cutting increased. 
Combining these new variables with those already selected, another 
set of screening runs was made from which the most promising equations 
were selected and the regression coefficients determined. Examination 
of these coefficients showed the signs on several of the transform- 
ations were not reasonable (including the site index transformations), 
and that a high degree of-multicollinearity existed. 
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Phase 4. Ridge regression was then used as a means of reducing the problems 
of multicollinearity. I hoped the minimization of multicollinearity 
might reverse the problem with sign on the site index transform, but 
the resulting analysis did not help. Because of the desire to have 
a model applicable over the range site indices naturally found in 
the habitat type, I decided to force a reasonable value of the site 
transform upon the model. This was done by defining a new dependent 
variable--the natural log of the quotient basal area growth divided 
by site index. By using this dependent variable, I assumed, therefore, 
that basal area growth would increase in a direct proportion to an 
increase in site index. 

Phase 5. A final set of screening runs was then made with the new dependent 
variable using the same methods described in the third phase of the 
analysis. From these screening runs, I selected promising models 
and made additional runs to determine their regression coefficients. 

Phase 6. Ridge regression was then used to selectively remove several of the 
time-since-last-cutting independent variables in order to minimize 
multicollinearity. Parameters to the final models were then deter- 
mined by ordinary, least squares regression. 


Log of Basal Area Growth Equations 


The finished log of basal area growth equations is in table 2. The two blackjack pine 
equations differ because of the presence or absence of the even-aged data when they were 
developed. The number of observations, the relative mean square residual (RMSQR),! and the 
coefficient of determination (R2 ) of each equation are found in table 3. 


Table 2.--Log of basal area growth equations 


Blackjack pine using all data: 


In(Basal Area Growth) = -8.51836897 + 1.16754330(1n(D)) - 4.00970143E-02(D) 
- 3.84298771E-03(LBAy) - 7.15483662E-03(MBA>) - 1.58234269E-02(UBA>) 
- 3.26097273E-01(A,1n(D)) + 8.80676713E-01(A3) + 1.0(1n(S)) 


Blackjack pine using uneven-aged data: 


In(Basal Area Growth) = -8,45357088 + 1.18165715(1n(D)) - 4.77068027E-02 (D) 
- 1.25420926E-03(LBAz) - 7.69300113E-03(MBA>) - 1.74056839E-02(UBA>) 
- 3.13247572E-01(A,1n(D)) + 9.10604169E-01(A3) + 1.0(in(S)) 


Yellow pine: 


In(Basal Area Growth) = -15.2464932 + 4.27656656(In(D)) - 1.83161626E-01 (D) 
- 7.27361567E-05(LBAz)2 - 9.11165626E-04(MBA,)2 - 2.41462106E-04 (UBA>)2 
- 1.05776062(A,) + 1.0(in(S)) 


where 
D = diameter class size 
S = Minor's (1964) site index 
MBA>y = basal area in the given diameter class plus the two adjoining 
larger and smaller diameter classes 
LBAy = total basal area below the smallest diameter class in MBA> 
= total basal area above the largest diameter class in MBA» 
A, = -0.244178 + 1.244178*EXP(-(1.176471 - 0.019607*TIME) 2) 
A3 = -0.000203 + 1.000203*EXP(-(1.428571 - 0.023809*TIME) ©) 
TIME = number of 5-year growth periods since last cutting 


1Relative mean square residual is the mean square residual for the model divided by the 
corrected variance of the dependent variable. See appendix C for more detail. 
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Table 3.--Statistics for log of basal area growth equations 


Equation Number of RMSQR R2 Modified Modified 
type observations RMSQR R2 

Blackjack pine 

with even-aged 

data S35 155 0.8001 0.2000 0.5155 0.4839 
Blackjack pine 

without even-aged 

data 27,268 .8223 STS) ~5425 .4576 
Yellow pine 4,558 . 9334 .0678 .8847 GS 


Because of the presence of within diameter class variation, the statistics quantifying the 
fit of the equations are unduly pessimistic. The equations explain between class variation 
and, therefore, the statistics were modified by eliminating the within diameter class variation 
component to better reflect the true fit of the equations. These "modified" statistics are 
also found) an table esr 


The results show that the blackjack pine equations developed with the even- and uneven- 
aged data explained over 48 percent of the between diameter class variation, the other black- 
jack pine equation explained over 45 percent, and the yellow pine equation explained only 
about 12 percent. For the blackjack pine, however, the equations do better at predicting 
managed growth than virgin growth. I obtained the statistics in table 4 by computing the 
residuals of each indicated data set around the final models. These statistics show that the 
blackjack pine equation developed with the even- and uneven-aged data, explained only about 19 
percent of the virgin uneven-aged variation, but it explained over 63 percent and almost 78 
percent of the variation for the managed uneven- and even-aged data sets, respectively. The 
other blackjack pine equations explained a little over 19 percent of the virgin and almost 64 
percent of the managed uneven-aged variation. Because the managed stand is of primary concern, 
these results are important. 


Table 4.--Modified coefficients of determination for various data sets by equation type 


Modified R2 


Equation Virgin Managed Managed 
type uneven-aged uneven-aged even-aged 
data data data 


Blackjack pine 


with even-aged data 0.1941 0.6334 0.7780 
Blackjack pine 

without even-aged data GNQZZ 6364 Sc 
Yellow pine .0440 . 0873 =e 


The poor results for yellow pine were not totally unexpected. These trees are overmature 
and frequently afflicted by various damaging agents that influence growth rates. In addition, 
the relatively small number of yellow pine in the stand often meant that between diameter 
class differences are actually between individual tree differences. Two mitigating factors, 
however, argued for retaining the equation. First, the model behaves as expected and does 
predict somewhat better than the mean. Second, the accurate prediction of yellow pine is not 
very important in managed stands because, under management, the presence of yellow pine usually 
would not be a favored condition. 


Check for Normality and Homogeneous Variance 


A check of the final equations revealed that the residuals were highly skewed and lepto- 
kurtic, which was the same result as first found for the individual plot equations (appendix 
C). Homogeneity of variance was checked by dividing the range of predicted log of basal area 
growth into intervals, computing a variance for each interval, and then comparing the variances 
across intervals for trends. Bartlett's chi-square test for homogeneity of variances (Snedecor 
and Cochran 1967) was not used because the test is highly sensitive to nonnormality. The 
visual checks for homogeneity indicated no consistent trends, and so weighting was judged 
unnecessary. 


Log Bias 


The usage of the log model to predict basal area growth introduces a problem because the 
value of interest for predictive purposes is basal area growth and not the log of basal area 


growth. The problem arises because EXP[in(Y)] is not a "mean-unbiased" estimator of Y; i.e.: 
aN 
EXP[E[in(Y)]] # E[Y] 
LIS - - - - 
Rather, EXP[1n(Y)] is a '"median-unbiased" estimator of Y (Bradu and Mundlak 1970). 


Because the residuals are not normally distributed, the log bias correction factors 
proposed by Bradu and Mundlak (1970), Oldham (1965), and Baskerville (1972) could not be 
used. Therefore, I proposed an alternative log bias correction factor that was added to the 
intercept term of the log of basal area growth model and was computed as the difference between 
the log of mean basal area growth minus the mean log of basal area growth (see appendix D). 


Models of Mean Residuals 


Because the residuals are not normally distributed, the mean deviation of the fast, 
moderate, and slow growers had to be computed empirically. I hypothesized that the range in 
residuals would increase as the time-since-last-cutting increased because cutting would, if 
properly applied, homogenize the stand by eliminating rough, cull, suppressed, and damaged 
trees. The residuals were, therefore, divided into time-since-last-cutting classes, and for 
each class the means of the upper, middle, and lower one-third residuals were computed. As 
expected, the means had a tendency to increase as time-since-last-cutting increased. Therefore, 
the time-since-last-cutting transforms (A,, Ap, and A3) were used to develop models that 
predict mean residuals for the upper and middle thirds as a function of time-since-last-cutting. 
The lower one-third mean model was then expressed as the negative of the sum of the other two 
models. This assumed that, for every value of time-since-last-cutting, the means of the three 
will sum to zero. The final residual models are in table 5. 


Table 5.--Models for the mean of upper, middle, and lower one-third of residuals about the 
indicated log of basal area growth model 


Blackjack pine with all data: 


Upper Mean = +0.56557130 + 0.12995917A, 
Middle Mean = +0.09120012 + 0.12946423A, 
Lower Mean = -0.65677142 - 0.25924340A, 


Blackjack pine with uneven-aged data: 
Upper Mean = +0.61870031 + 0.08118138A3, 


Middle Mean = +0.13230685 + 0.08389430A3 
Lower Mean = -0.75100716 - 0.16507568A3 


+ + 


Yellow pine: 


Upper Mean = +0.70338067 
Middle Mean = +0.24266226 
Lower Mean = -0.94604293 


as 


0.7280557A3 
0.59168455A, 
1.31974025A3 


+ 


where 
A, = -0.244178 + 1.244178*EXP(-(1.176471 - 0.019607*TIME) 3) 
Ay = -0.095336 + 1.095336*EXP(-(1.250000 - 0.020833*TIME) +) 
A3 = -0.000203 + 1.000203*EXP(-(1.428571 - 0.023809*TIME) ©) 


Prediction of Basal Area Growth and Diameter Growth 


Basal area growth of the one-third fastest (BAG,), moderate (BAG»), and slowest (BAG3) 
growing trees is predicted by: 


BAG, = EXP[R, + K + In(BAG)], i = 1, 2, or 3 


where: 
BAG . = predicted basal area growth for the ith third of the residuals 
R. = predicted mean log of basal area growth residual for the ith 
third of the residuals 
K *= log bias correction) factor 
In(BAG) = final, average log of basal growth model 


Diameter growth of each third can then be computed by the relationship 


DG, = yo (DBA + BAG.) - q22DBA 
i TT i TT 


where: 
DG, = predicted diameter growth for the ith third of the residuals 
D = diameter class size 


DBA = average basal area of the diameter class 


T 


576 


(D2 + 0.1D + 0.085) 
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The equation for the average basal area of a diameter class assumes that the trees in a 
diameter class are uniformly distributed over the 1/10-inch subclasses of the class. As 
previously shown, this assumption could not be disproven and therefore would appear correct. 


An Initial Evaluation of the Log of Basal Area Growth Equations 


An unanswered question is whether basal area growth for even- and uneven-aged stands is 
logically interrelated. As mentioned previously, the problems with the data used in this 
study and the lack of normality make a statistical answer to this question intractable. I 
thought, however, that an initial comparison between the two blackjack pine equations could 
be made in such a fashion as to shed light on the performance of each model. A second compari- 
son and a selection of a final blackjack pine model were reserved until validation. 


To make this initial comparison, the data used to develop the equations were divided into 
four sets: virgin uneven-aged data, managed uneven-aged data, managed even- and uneven-aged 
data, and managed even-aged data. The residuals about each final log of basal area growth 
equation were then computed separately for each of the data sets originally used to develop 
the particular equation. From this information, a mean squared error (MSE), average residual, 
and percent average residual (expressed as percent of the mean log of basal area growth) were 
computed for each data set. The results are in tables 6 through 9. 


Table 6.--Comparison of actual log of basal area growth data from virgin, uneven-aged plots 
to log of basal area growth predicted by various equations developed using data 
from those plots indicated 


Plots used Number of Mean squared Average Percent 
Vigor class in equation observations error residual average 

deve lopment residual 
Blackjack pine All 8,197 AS72575 +0.000013 -0.000462 
Blackjack pine Uneven-aged 8 h97 7 5900 + .000028 - .000994 
Blackjack pine Virgin, 

uneven-aged 85197 1.534644 +0 0 
Yellow pine Uneven-aged 3,762 6.310643 + .000261 - .000076 
Yellow pine Virgin, 

uneven-aged 3,02 6.273243 +0 0 


Table 7.--Comparison of actual log of basal area growth data from managed, uneven-aged plots 
to log of basal area growth predicted by various equations developed using data from 
those plots indicated 


Plots used Number of Mean squared Average Percent 
Vigor class in equation observations error residual average 

development residual 
Blackjack pine All 19,069 0.940228 -0.007303 +0.257221 
Blackjack pine Uneven-aged 19,069 - 938867 - .000012 + .000423 
Blackjack pine Managed, 

uneven-aged 19,069 935147 0 0 
Yellow pine Uneven-aged 794 1.840926 - .001237 + .058064 
Yellow pine Managed, 

uneven-aged 794 1.800258 0 0 
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Table 8.--Comparison of actual log of basal area growth data from managed, even-aged plots to 
log of basal area growth predicted by various equations developed using data from 
those plots indicated 


Plots used Number of Mean squared Average Percent 
Vigor class in equation observations ERO residual average 
development residual 
Blackjack pine All 5,886 0.224909 +0.023643 -0.814989 
Blackjack pine Managed, 
even-aged 5,886 hoi Sil 0 0 


Table 9.--Comparison of actual log of basal area growth data from managed, even- and uneven- 
aged plots to log of basal area growth predicted by various equations developed 
using data from those plots indicated 


Plots used Number of Mean squared Averege Percent 
Vigor class in equation observations error residval average 
development residual 
Blackjack pine All 24,956 OL 7713355 -0.000004 +0. 000140 
Blackjack pine All 
managed 24,956 . 767860 0 0 


Also in those tables are the MSE's, average residuals, and percent average residuals for 
selected equations developed in phases 2 and 3 of the analysis using just the data from the 
given set. The equations chosen for comparison were those minimizing RMSQR for each data set. 
In choosing these equations, no concern was given to whether they were reasonable in behavior, 
and most of them were not. These equations do, however, represent the best that could be done 
for the given data set and, therefore, provide a benchmark for how good the final equations 
are at predicting growth. 


For blackjack pine, the results for the virgin uneven-aged data (table 6) indicate that 
the two equations are quite similar at predicting growth. For the managed, uneven-aged data 
(table 7), the equation developed with the even-aged data shows a larger (but still small) 
average residual, but the MSE's are almost the same. This positive percent residual is offset 
by a negative percent residual for the even-aged data (table 8). When the managed even- and 
uneven-aged data are put together (table 9), the average residual becomes insignificant. The 
conclusion, therefore, is that the two blackjack equations give about the same results as far 
as MSE's and average residuals are concerned. 


For yellow pine, the results indicate not much predictive capability was lost when the 
two data sets (virgin and managed uneven-aged data sets) were combined. 


Mortality 


Mortality is often separated into a catastrophic (irregular) or a noncatastrophic (regular) 
form (Lee 1971; Stage 1973; Monserud 1976). As its name implies, catastrophic mortality is 
the result of usually unpredictable, massive disturbances such as fire, hurricanes, tornadoes, 
and insect or disease epidemics. Noncatastrophic mortality accompanies normal stand develop- 
ment. The models developed in this simulator predict only noncatastrophic mortality. 
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The following is a summary of the procedures used to obtain these models. Additional 
information is in appendix E. 


Prior silvicultural and mensurational research in southwestern ponderosa pine showed that 
noncatastrophic mortality in uneven-aged stands is often correlated to diameter size, severity 
of cutting, and tree vigor as indicated by bark color (in other words, whether the tree is a 
blackjack or yellow pine). For even-aged stands, site and total stand basal area are signifi- 
cant in predicting total stand mortality. 


I selected a nonlinear logistic function for modeling the proportion of trees in a dia- 
meter class succumbing to mortality. One reason for this is that the dependent variable of 
the logistic can be properly limited to a value between zero and one. Second, characterizing 
a dichotomous dependent variable such as mortality with the logistic function appears to 
produce a model with improved statistical properties when compared to ordinary, linear least 
squares models (Hamilton 1974; Hamilton and Edwards 1976). 


To model the proportion of the trees in a diameter class dying in a growth period, the 
dependent variable was set to zero if the tree survived to the end of the growth period, and it 
was set to one if the tree died. The independent variables used were transformations of 
diameter class or stand attributes and incorporated not only the prior factors found to be 
correlated to mortality, but also factors hypothesized to be correlated. Included in this 
group were time-since-last-cutting and basal area growth predicted from the three completed 
models. Appropriate transformations were formed based upon both prior silvicultural and 
mensurational research findings, and an examination of the particular behavior of the logistic 
function. 


The development of finished mortality models required the use of nonlinear regression and 
proceeded in three phases: 


Phase 1.--For each vigor class, the data were separated into three sets: virgin uneven- 
aged, managed uneven-aged, and managed even-aged. Preliminary nonlinear regression runs, 
using the logistic function, were made to select those variables most highly correlated to 
mortality. The two best variables proved to be diameter class size squared and predicted 
basal area growth. 


Phase 2.--The separated data sets were collapsed into two sets (uneven-aged, and even- 
and uneven-aged) by modeling the change in model parameters between the separated data sets as 
functions of time-since-last-cutting. An examination of a chi-square "goodness-of-fit" test 
for the best models indicated that the uneven-aged yellow pine equation (table 10) and the 
uneven-aged blackjack pine equations (tables 11 and 12) were satisfactory. The equations for 
even- and uneven-aged blackjack pine were not satisfactory. 


Phase 3.--To resolve the problems with the even- and uneven-aged blackjack pine data set, 
another set of nonlinear regression runs was made. While these models were an improvement 
over the previous models, the chi-square values still indicated a significant difference 
between predicted and actual mortality (tables 13 and 14). An examination of the causes for 
mortality revealed that the even-aged data had incurred considerable snowbreak loss. Because 
I could not decide whether snowbreak should be treated as a catastrophic or noncatastrophic 
loss, these models were retained until the validation phase of the study could more closely 
examine the predictive capabilities of the equations. 


13 


Table 10.--Chi-square test across diameter classes for best mortality equation for yellow pine 


Diameter Number of trees Actual Predicted Chi-square 


class in class mortality mortality value 
4 - 11 83 2 Wee 0.0 
12 - 13 12/9 5) Pehle 0 
14 - 15 145 3 2.79 0 
N69= 17 229 5 ANI, ae 
18 - 19 525 5 5.42 0 
20 - 21 469 4 USL 1S) 
22 = 123, 585 6 8.93 1.0 
24 - 25 630 15 E56 Spill 
26 27, 557 4 8.68 Das) 
28 - 29 535 13 9.14 1.6 
20) soil 370 6 6.83 gil 
BD = 53 260 10 SenS Soe 
34 - 35 154 5 3.84 we. 
SO. = 9 Sif) 105 2 3.03 .4 
38 = 39 47 2 59 gal 
40+ 30 0 1.25 1S 


|: 


— 


Chi-square statistic = 15. 


Table 11.--Chi-square test across diameter classes for best mortality equation for uneven-aged 
blackjack pine using basal area growth equation for even- and uneven-aged blackjack pine 


Diameter Number of trees Actual Predicted Chi-square 
class in class mortality mortality value 
4 4,416 11 OFA 3.4 
5 4,066 18 WSF, 0 
6 $4522 18 AR heal 
7 2,629 17 TASS 2.8 
8 20211 8 8.81 gl 
9 Soy 7 6.82 0 
10 AOS) ) 5259 0 
11 1,004 6 4.50 35 
12 893 3 4.01 52) 
13 737 2 3.52 aes) 
14 686 1 Sus INGE 
15 666 6 Sl9 Dos 
16 685 4 51150 oll 
7, 649 2 Sy IY .4 
18 610 1 3.06 1.4 
19 473 0 2x55) Zin 
20 - 21 908 6 4.83 a5 
22. =) (25) 622 9 3.41 92 
242 = 225 353 1 195 oe) 
26+ 261 0 1.94 139 


oO 


Chi-square statistic = 28. 
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Table 12.--Chi-square test across diameter classes for best mortality equation for uneven-aged 
blackjack pine using basal area growth equation for uneven-aged blackjack pine only 


Diameter Number of trees Actual Predicted Chi-square 


class in class mortality mortality value 
4 4,416 11 19.10 3.4 
5 4,066 18 17.36 0 
6 53522 18 14.11 ue 
ia 2629 17 11.32 2.8 
8 2,021 8 8.81 oll 
a 1 So 7 6.81 0 
10 1,203 5 5.59 0 
ll 1,004 6 4.50 <5 
12 893 3 4.02 +3 
13 737 2 5/352 25 
14 686 1 See) 1.5 
15 666 6 3.20 25 
16 685 4 3 Sh zl 
17 649 2 Bo Ks) -4 
18 610 1 3-07 1.4 
19 473 0 Des4 2235 
20 - 21 908 6 4.84 ac) 
22 - 23 633 e) 3.41 7 
24 = 25 353 1 WG. 35 
26+ 261 0 195 19 


Chi-square statistic = 28.9 


Table 13.--Chi-square test across diameter classes for best mortality equation for even- and 
uneven-aged blackjack pine using basal area growth equations for even- and uneven- 
aged blackjack pine 


Diameter Number of trees Actual Predicted Chi-square 
class in class mortality mortality value 
4 5,790 126 83.78 2125 
5 5,635 47 70.05 7.6 
6 4,681 28 44.95 6.4 
7 3,549 22 26.04 -6 
8 2,540 8 14.34 2E8 
s) 1,792 7 8.56 S 
10 1,261 6 5.68 0 
1l 10nd 6 4.23 a7 / 
12 893 3 3.54 oil 
13 DSit, 2 2-95 2 
14 686 1 25S =9 
15 666 6 2.61 4.4 
16 685 4 2.63 7 
LF, 649 2 2252 al 
18 610 1 2.38 .8 
19 473 0 1.96 2.0 
20 - 21 908 6 4.07 -9 
22.25 622 9 3.06 TS 
24 - 25 353 i PANES -6 
26+ 261 0 2.50 2S 


Chi-square statistic = 64.53 


a aa 
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Table 14.--Chi-square test across diameter classes for best mortality equation for even- and 
uneven-aged blackjack pine using basal area growth equations for uneven-aged black- 
jack pine only 


Diameter... +») ~Number of trees) Actually Salm) mPredilctediy mn = aNNChi=squancum 
class iim cllasis mortality mortality value 
4 520, 126 83.56 Pla 
5 5,635 47 70.76 8.0 
6 4,681 28 44.83 Ons 
7 3,549 22 29.458 5 
8 2,540 8 13.96 7S 5 
9 1,792 7 8.45 2 
10 12611 6 Sal 0 
TAL Oa: 6 4.30 wi 
12 893 5 SoZ il 
13 U5 2 ZO ae 
14 686 1 Poy 0 
LS 666 6 Dif dl 4.0 
16 685 4 Pot .6 
Ly 649 2 2.63 52 
18 610 1 2.49 >) 
Me) 473 0) 2305 Doll 
20 - 21 908 6 AyaS .8 
22a" 25 622 9 S705 SS 
24 - 25 555 1 2.04 5.) 
26+ 261 0 2.34 De® 


Ww 


Chi-square statistic = 64. 


The finished mortality equations are in table 15. Graphs for some of these equations, using 
representative values of thelandependent variables, areian filgunes 15 25 vand) 5. Duew covthe 
Similarity of the two basal area growth equations for blackjack pine, only one of the two 
mortality equations for uneven-aged blackjack pine and for even- and uneven-aged blackjack 
pine were plotted. In all of the equations, the mortality rate increases as D and MBA») 
increases and as BAG and TIME decreases. Caution should be taken when interpreting these 
graphs. Because of the interrelationships among D, BAG, and MBA,, the actual mortality rate in 
a stand would not be represented by any single curve shown. 


Conversion from Blackjack Pine to Yellow Pine 


The classification of ponderosa pine as either blackjack pine or yellow pine is based on 
bark color. Young, vigorous growing trees have dark-colored bark while the mature or over- 
mature, slow growing trees develop a yellow-colored bark (Harlow and Harrar 1958). As this 
study has shown, the trees in these two ''vigor'' classes display different growth and mortality 
rates and, as a result, an improvement in modeling growth and mortality can be expected if the 
two vigor classes are treated separately. By treating them separately, however, it then 
becomes necessary to develop a model predicting the conversion from blackjack pine to yellow 
pine. 
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Table 15.--Finished mortality equations 


PM = (1.0 +e) 


where 


PM = proportion of trees in a diameter class dying in the next 5-year period 


For uneven-aged yellow pine: 


X = -3.612474 - 1.763910*YPBAG - 9.281652*A,*YPBAG + 5 .597303E-04*D2 


For uneven-aged blackjack pine: 


X -5.372649 + 1.270862E-03*D2 - 9.599418*A, *BJBAG] or 


X 


-5.372991 + 1.271523E-03*D2 - 9.388071*A, *BJBAG2 


For even- and uneven-aged blackjack pine: 


-2 
-4.670214 - 18.54357*BJBAG1 - 7.948837E-03*TIME + 1.870818* (TIME) 
+ 1.370731E-02*MBA> + 4. 382886E-03*D? or 


X 


-2 
-4.725930 - 17.46734*BJBAG2 - 7.256918E-03*TIME + 1.883188* (TIME) 
+ 1.423927E-02*MBA> + 4.168657E-03*D2 


Past 
iT) 


D = diameter class size 
YPBAG = predicted yellow pine basal area growth 
BJBAG] = predicted even- and uneven-aged blackjack pine basal area growth 
BJBAG2 = predicted uneven-aged blackjack pine basal area growth 
TIME = number of 5-year periods since last cutting 
A,, Ao, A3 = sigmoidal transforms of TIME. See table 5 for equations. 


MBAz = basal area per acre in the given diameter class plus the two adjoining larger 
and smaller diameter classes 
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Figure 1.--Predicted 5-year mortality rate for yellow pine 
diameter class size and basal area growth (YPBAG). 
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Figure 2.--Predicted 5-year mortality rate for uneven-aged 
blackjack pine by diameter class size and uneven-aged basal 
area growth (BJBAG2). 
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Figure 3.--Predicted 5-year mortality rate for even- and uneven-aged blackjack pine by diameter 
size class, MBA» class, and even- and uneven-aged blackjack pine basal area growth (BJBAG1). 


During data collection, all trees were initially categorized by vigor class, and then 
those trees not cut on a plot were reclassified every 20 years (that is, at the start of 
every fifth growth period). An examination of the data available for modeling conversion 
indicated the following problems: 


1. While growth, mortality, and ingrowth could be predicted every growth period, conver- 
sion from blackjack to yellow pine could only be predicted after every fourth growth period. 
This might cause the number of trees in each vigor class to change abruptly after every fourth 
growth period. 


2. Cutting of blackjack trees between the reclassification times negates the utility of 
the plot for predicting conversion because of the uncertainty of whether those trees cut would 
have converted or not. The plot and classification period combination without cutting are the 
first classification period of plot 62 and both periods of plot 61. 


3. The conversion from a blackjack pine to a yellow pine is a slow but continual process, 
and therefore, the discrete classification of a tree as either blackjack or yellow pine is 
extremely subjective. This subjectivity can become a problem on long-term studies because the 
classification of vigor classes is done by different people with individualized critera for 
classification. An examination of the conversion rates for the first 20 years on plots 61 and 
62, compared to the conversion rate on plot 61 in the second 20 years, shows that the rates 
differ considerably between the two classification periods. While it cannot be proved, it is 
likely this difference is due in part to the different individuals performing the classifica- 
tions. 


Based on these problems, I decided to use only the data from plot 61 (both classification 
periods) to develop the conversion rate model. The data for the first classification period 
on plot 62 were eliminated since including that data would unduly weight the resulting model 
towards the low conversion rate of the first classification period. By using just the data 
from plot 61, an approximately equal weight is given to both periods. 
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The usage of data from the second classification period on plot 61 does pose an addi- 
tional problem. During the second classification period, only trees 7.6 inches and larger 
were measured. All stand density variables, therefore, must be expressed in trees 7.6 inches 
and larger. In addition, the usage of predicted basal area growth is precluded because of the 
equation's dependency upon knowing stand basal area of all trees 3.6 inches and larger. 


The following is a summary of the process used to develop the final conversion model. 
Further details are in appendix F. 


I anticipated the factors influencing the conversion rate would be the same as those that 
might influence tree vigor, as expressed by growth. Previously it was shown these factors 
include measures of productivity, diameter class size, and stand density. Keeping in mind 
both these factors and the restriction that only stand attributes for trees 7.6 inches and 
larger could be used, a number of variables were defined for use as independent variables. 


As with mortality, I decided to express the rate as a proportion of the trees in a diameter 
class converting from blackjack to yellow pine and to model this rate using the nonlinear 
logistic function. Therefore, I assigned the dependent variable a value of one if the black- 
jack tree converted to yellow pine, and a value of zero if it did not. 


Using these variables, a number of nonlinear regression runs were made, and chi-square 
"goodness-of-fit"' values were examined. Predictions from all models displayed significant 
differences from the actual values. Examination of the models and data revealed that the 
greatest misfit was due to data from a single subplot. Therefore, I decided to strengthen 
the data base by adding to it the data from those subplots originally eliminated because of 
the paved highway. A second set of nonlinear regression runs was made, but again the resulting 
models displayed poor fit in the same diameter class. 


Because tree vigor is such a subjective classification, the apparently atypical subplot 
was eliminated and new equations determined. The resulting model's predictions were not 
Significantly different from the actual values (see table 16). To test this model further, 
chi-square values were also computed for the data reserved for validation. The result indicated 
a significant difference between the model predictions and the actual data, and the cause was 
traced to another subplot with an unusual conversion rate in the same diameter class. I could 
not determine why the two subplots exhibited such high conversion in the same diameter class, 
but the elimination of the second subplot from the data did produce an insignificant chi- 
square value for the validation data. 

The final conversion model is: 


eyes 
PrCONe= (Gl302- E>) 


where 
2 
X = -3.777269 + 0.5012605D - 6.597465E-03D - 0.0982176S 


PrCON = proportion of trees in a diameter class converting from blackjack pine to yellow 
pine 


D = diameter class size 


S = Minor's site index 

A plot of this equation is in figure 4. To use this model, each independent variable is 
assigned the value existing for the diameter class at the start of the growth period prior to 
reclassification. Therefore, at the start of every fourth growth period, a prediction is made 
as to how many trees will convert by the start of the fifth growth period. 


PAL 


Table 16.--Chi-square test across diameter classes for final conversion model, developed using 
the expanded data set, without data from subplot 13 


Diameter Number of trees Actual Predicted Chi-square 
class in class conversion conversion value 
A= 9 1,486 0.00 53 1S 
10 - 11 317 0 .94 59) 
12S, 228 1.00 59) al 
14 103 4.00 1.06 Sod 
1S} 99 1.00 56 vl 
16 91 6.00 Wai HORS 
17 91 1.00 2.43 sfs) 
18 93 4.00 2.81 xo 
19 100 4.00 4.03 0 
20 87 2.00 4.70 AS 
21 79 5.00 4.87 0 
22 74 4.00 5629 Be 
BS 64 8.00 oIk7/ 5) 
24 47 1.00 Syaz 3.4 
25 42 6.00 5.69 0 
26 52 4.00 Seles ae 
27 24 5.00 4.13 ot 
28 15 4.00 5.08 el 
29+ 19 5.00 4.32 bal 


Chi-square statistic = 28.7 
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Proportion of Trees in a Diameter Class Converting in Twenty Years 


Figure 4.--Predicted conversion 
from blackjack pine to yellow 
pine by diameter class size 
and site index (S). Diameter Class Size in Inches 


4 8 12 16 20 24 


Recruitment 


The requirements for adequate natural regeneration in the southwestern ponderosa pine are 
the presence of large, vigorous trees well dispersed over the area, a good seed crop, adequate 
rainfall in two consecutive years, and low level of competitive vegetation (see appendix G for 
details). This combination of conditions occurs at erratic invervals and, as a result, ade- 
quate natural regeneration is erratic. For example, the last excellent regeneration period 
occurred in 1918-1919 (Pearson 1950; Meagher 1950; and Schubert 1974). Cooper (1960) has 
concluded that erratic natural regeneration is a major reason for the "unbalanced" size struc- 
ture of the ponderosa pine forests of the Southwest. 

Ideally, a natural regeneration model should take all of these factors into consideration 
as components of it (such as the one described by Leak and Graber 1976). Unfortunately, the 
data available in this study preclude the development of a regeneration model. This was 
unfortunate because the greatest distinction between even- and uneven-aged stand development 
is the method of regeneration. The even-aged stand usually starts as a surge of regeneration 
over a relatively short period of time (20 years or less). Next is a stage where little or no 
regeneration takes place. Finally, some regeneration may again occur (depending upon the 
species and site characteristics) as the stand becomes mature (or overmature) and develops a 
more open structure. On the other hand, the uneven-aged stand is usually characterized by a 
more or less continuous flow of regeneration into the stand. 
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Because of the data restrictions, the only “recruitment! model possible atrthasmermesas 
an ingrowth model. In this study, ingrowth is defined as the number of new trees growing into 
the 4- or 5-inch diameter classes during a growth period. Ingrowth, therefore, includes not 
only those factors influencing the amount of established seedlings, but also includes those 
factors influencing the growth and mortality rate of seedlings as they advance to ingrowth 
size. Conditions providing for a large, established seedling crop and allowing favorable 
seedling growth and low mortality will ultimately result in a large ingrowth rate at some 
point ine the- future: 


For uneven-aged stands, ingrowth can be expected at any time, but in even-aged stands, 
ingrowth would only be expected during the early, and possibly the late, stand development 
phase. Unfortunately, the Taylor Woods even-aged data represent only stands in relatively 
early stand development. As a result, any ingrowth model developed using just the Taylor 
Woods data would be applicable only to stands in the same stand development phase. Such a 
restriction would give the model little utility for predicting even-aged stand development 
through all phases of development. Therefore, it was decided that the ingrowth model would be 
developed only for uneven-aged stands. An alternative method for even-aged stands by which 
use of the basal area and mortality models are extended down into the l1-, 2-, and 3-inch 
diameter classes is proposed and tested in the validation phase. 


Two equations were used to model uneven-aged ingrowth. The first predicts the total 
number of ingrowth trees. The second predicts the proportion of the total ingrowth trees that 
will grow through the 4-inch diameter class and into the 5-inch diameter class during the 
growth period. 


TOTAL INGROWTH 


After examining prior silvicultural and mensurational findings, I hypothesized that the 
factors possibly influencing the rate of ingrowth for uneven-aged stands included potential 
of the stand to produce cones, level and structure of competition, site index, mean stand 
diameter, and quadratic mean stand diameter. Based on these factors, appropriate independent 
variables were defined. The following is a summary of the process used to develop the total 
ingrowth model. Further details are in appendix G. 


Two factors influenced the development of a strategy for modeling total. ingrowth. First, 
while a nonlinear total ingrowth model is most likely, it was not known whether the residuals 
about the model are multiplicative (implying that taking the log of the model could linearize 
the data) or additive. Second, the form and specification of the applicable independent 
variables in the final ingrowth model had not been rigorously defined. Therefore, a good deal 
of screening was necessary before the final model could be determined. Given these conditions, 
the following strategy was used to develop a total ingrowth model: 


Phase 1.--The total ingrowth data were separated into two sets: virgin uneven-aged and 
managed uneven-aged. For each set, a preliminary all-combinations screening run using linear 
regression was made on the variables linearized by applying logarithms, and selected models 
were chosen for further analysis. 


Phase 2.--A second screening run combined the two basic data sets through use of the 
time-since-last-cutting transforms. Because the sign on site index was not considered reason- 
able, a new dependent variable was formed by dividing total ingrowth by site index and taking 
the log of the quotient. A new set of screening runs was then made and the final independent 
variables were chosen. 


Phase 3.--An examination of the residuals of the preceding model indicated that the error 
structure of total ingrowth divided by site index might not be best represented by a log 
model. Therefore, the antilog of the model was fitted with a linear, least squares slope 
correction. An examination of these residuals revealed that the error structure of total 
ingrowth divided by site index was better represented by a weighted, nonlinear model with an 
additive error. 
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Phase 4.--Weights were then developed and used in a weighted, nonlinear regression 
program to determine the following, final weighted nonlinear model: 


di, di, 
Total ingrowth = gS (BACH )O (1. O+BACH,)2e003 + dsA;BA *) 


where 
dg = 0.48630258 d3 = -3.97721395E-06 
d,; = 0.65949418 d, = 2.7667388 
d3 = -0.49280582 ds = -1.33299951E-06 
where 


S = Minor"s (1964) site index 

BA = total stand basal area per acre in square feet 

BACL, = basal area per acre in the 4- through 6-inch diameter classes 
BACL» = basal area per acre in the 7 -inch diameter classes. 


The weighted, mean square error for this model was 0.78244674. The coefficient of determina- 
tion, R2, was 0.8285 for 201 observations. See figure 5 for a graph of this function. 
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Figure 5.--Predicted total number of ingrowth trees by site index classes (S), 
total basal area classes (BA), and the basal area in the 4-, 5-, and 6-inch 
diameter classes. 
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THROUGH-GROWTH 


Through-growth is defined as the proportion of the total ingrowth that will grow through 
the 4-inch diameter class and into the 5-inch class, in the same growth period. It is assumed 
that no trees will grow through to the 6-inch class. 


As with the mortality and conversion models, I modeled through-growth using the nonlinear 
logistic function. Independent variables tried in this weighted, nonlinear analysis included 
predicted basal area growth in the 4-inch diameter class, predicted total ingrowth, and time- 
since-last-cutting. The finished models were selected after fitting and closely examining a 
number of alternative models. The regression coefficients for these models are in table 17 
and plots of the models in figures 6 and 7. 


Table 17.--Final through-growth equations 


Through-growth (1.0 +e 7) 
where 


Through-growth = proportion of the total number of ingrowth trees that through grow into 
the 5-inch diameter class 


-2.636336 + 28.73438*BJBAG1-4" - 0.03419626*A, *TINGRO 
Xx = or 
-2.490266 + 24.95862*BJBAG2-4" - 0.03440381*A, *TINGRO 


BJBAG1-4'"' = predicted basal area growth of the 4-inch diameter class using the even- and 
uneven-aged blackjack pine equation 


BJBAG2-4'' = predicted basal area growth of the 4-inch diameter class using the uneven- 
aged blackjack pine equation 


TINGRO = predicted total number of ingrowth trees 
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Figure 6.--Predicted through-growth by 
BJBAG1-4"" and total ingrowth (TI). 
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Figure 7.--Predicted through-growth by 


BJBAG2-4"" and total ingrowth (TI). 


Height and Cubic-foot Volume Equations 


The lack of height-growth data necessitated developing alternative height equations for 
assessing product potential as predicted by cubic-foot volume equations. Two separate 
equations were developed for even-aged and for uneven-aged stands. 


EVEN-AGED HEIGHT EQUATION 


Minor (1964) developed an equation for even-aged blackjack stands in northern Arizona 
that predicts heights of dominant trees as a function of site index and stand age. Plots of 
average diameter class heights over diameter class size for the Taylor Woods height data 
indicated that average diameter class height was below the height predicted by Minor's (1964) 
equation for the smallest diameter classes and increased, asymptotically, to Minor's (1964) 
height prediction in the largest diameter classes. Using the procedure outlined in appendix 
H, the following equation was developed for use in even-aged, blackjack stands: 


H, = MH - 0.88015037(MH-4.5)(1.0 - ee 
where 
H. = average height of the Dth diameter class 
= diameter class size 
MH = maximum height as predicted by Minor's (1964) equation 
=Si A003 GVA = 10)) 10 1559(S)' (VA = 10) 
S = Minor's (1964) site index 


A = breast height stand age, 20<A<140 


DM = maximum diameter class size that the stand has achieved in its development. 


UNEVEN-AGED HEIGHT EQUATION 


Because stand age is meaningless in uneven-aged stands, the foregoing even-aged equation 
cannot be applied to uneven-aged stands. As an alternative, an equation that predicts average 
diameter class height as a function of site index and diameter class size was developed using 
the procedure described in appendix H. The equation is: 


(D+35)~ 2 
= Joy ar Jax) o 5°72 c 2O 


22) 
\ 


where 
Hp = average height of the Dth diameter class 
D = diameter class site 
S = Minor's (1964) site index 
bg = 4.5 
by = 13.178649 
by = 0.71631005 


b3 = -4221.6528 
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This equation is for both blackjack and yellow pine. Insufficient data precluded separate 
equations. 


TOTAL STEM CUBIC-FOOT VOLUME EQUATIONS 


Existing equations or data for developing new equations to predict total stem cubic foot 
volume in l-inch diameter classes were not available for northern Arizona. As an alternative, 
I used tree volume equations developed by Hann and Bare (1978) for yellow and blackjack pine 
on the Coconino National Forest of northern Arizona. While application of an individual tree 
volume equation to diameter class attributes could produce biased volume estimates, I felt the 
smallness of l-inch classes would minimize the potential problem. 


Structure of Simulator 


The four models--upgrowth, mortality, conversion from blackjack pine to yellow pine, and 
ingrowth--interact in the following fashion: 


1. Present number of trees in each diameter and vigor class provide the starting point. 


2. Mortality is computed and subtracted from each class to give the number of trees 
expected to survive to the end of the growth period. 


3. Using diameter class growth of the survivors and the appropriate within diameter 
class distribution, the number of trees moving to various diameter classes is determined. 


4. At the start of each fourth growth period, the blackjack-to-yellow pine conversion 
rates for survivors are calculated; otherwise, the conversion rates are set to zero. The 
conversion rates and upgrowth information then are used to allocate the survivors to the 
appropriate "future" diameter and vigor classes. 


5. Ingrowth for the next 5-year period is computed and added to the 4- and 5-inch classes. 


6. Removals due to cutting are made from each diameter class. Cutting, therefore, 
presumably occurs at the end of the growth period. This step completes the "future'' diameter 
distribution. 


7. This distribution then becomes’ the present: diameter idistribution, and ‘the cycle 
starts again. A system interaction chart is found in figure 8. 


VALIDATION OF SIMULATOR 


The problem of validating a simulator has troubled modelers for years. Shannon (1975) 
identified and described the philosophies of three approaches. "Rationalism holds that a 
model is simply a system of logical deductions from a set of premises, which may or may not be 
open to empirical verification or appeal to objective experience" (page 212). ''Empiricism 
refuses to admit any premises or assumptions that cannot be verified independently by experi- 
ment or analysis of empirical data'' (page 214). Finally, absolute pragmatism holds that a 
model is designed to meet a need, and if the model succeeds in meeting the need, then it has 
been validated. 
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Figure 8.--System interaction chart. 


People interested in validation have come to realize that none of the three philosophies 
fully meets validation needs individually or is applicable to all situations. As a result, a 
utilitarian approach has evolved. Shannon defines the three steps of this approach: 


The first stage is to seek face validity of the internal structure of the model 
based upon a priori knowledge, past research, and existing theory (page 215). 

The second stage is also concerned with the validation of the internal structure 
of the model, and consists in empirically testing, whenever possible, the hypothesis 
used (page 215). 

The third stage attempts vigorously to verify the model's ability to predict the 
behavior of the real world system (page 216). 


Many applications to forestry require the simulator to be both a prediction (or "prognosis"') 
tool and a technique for examining the dynamics of the system. This requires the simulator to 
be accurate not only in a behavioral sense but also quantitatively. Many of the system com- 
ponents are slow to react and, coupled with the general long-range aspects of forestry, require 
long simulation runs to fully monitor true system behavior. With these special characteristics 
in mind, I slightly modified and expanded upon the preceding utilitarian philosophy to produce 
the following tentative "rules of validation." 


1. The simulator components should have a sound theoretical basis that uses as much 
experimental evidence and expert knowledge as possible. This would include identifying the 
Significant system components, the factors that can influence these components, and applicable 
forms and techniques for modeling them. 
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2. Correct statistical and modeling techniques should be used to develop estimates of 
the parameters of the basic components. 


3. Where possible, the estimated parameters and model forms of these components should 
be validated with a data source independent from that first used to estimate the parameters. 


4, Like its components, the validation of the simulator should also be made on a data 
set independent of those used to derive the simulator. This is necessary only if the simulator 
is to be applied to studies other than those used in its development (which is the case most 
of the time in forestry applications). This or these independent data set(s) should also have 
a form compatible with the underlying structure of the simulator and should cover as wide a 
range of conditions as possible. 


5. Validation should occur on the basic attributes that the simulator produces. For 
example, if the simulator is an individual tree/distance dependent one, then validation should 
occur on the tree level attributes rather than on stand level attributes. Aggregation tends 
£0) cover up inconsistencies: 


6. Validation should occur over long enough time to allow valid tests (or comparisons) 
of long- as well as short-term responses and of subtle changes that take a long time to 
manifest themselves. 


7. Where possible, some type of ''statistical" test should be used to aid in quantifying 
validity. 


Applying the Rules 


The first three rules are normally applied during the model building process, while the 
last four rules are applicable to the final simulator. To incorporate the first two rules in 
this study, the following steps were taken in the model building process: (1) an exhaustive 
examination of the literature to identify major components of the simulator, the factors 
influencing these components, and applicable model forms; (2) an insistence that the various 
equations used to model these components meet the expected model forms and behave reasonably 
across the expected input data range; (3) identification of the correct statistical tools to 
model the components; and (4) where possible, a test of the assumption of both the models and 
Statistical tools: 


I decided, however, to develop a validation process that combined rule 3 with the last 
four rules so I could evaluate each component by how it influences the prediction of future 
diameter distributions, and by the necessity of saving time and conserving monetary resources. 
To combine rule 3 with the last four rules, it was necessary to make four different validation 
runs. The first run used actual values for all components except upgrowth;. that is, only 
upgrowth was predicted. The second run used predicted upgrowth and mortality. The third run 
used predicted upgrowth, mortality, and conversion. The final run predicted all components 
(upgrowth, mortality, ingrowth, and conversion). In this way, the effect of adding each 
predicted component to the simulator could be evaluated by how the simulator predicts the 
basic attribute of interest, the future diameter distributions (rule 5). 


Because two blackjack pine equations were developed (one with the even-aged data and one 
without), both equations were tested on the first set of validation runs. The equation that 
performed best was the final blackjack pine growth equation. A similar approach was used on 
the second set of runs to evaluate various mortality models. 


Validation runs were made on the data reserved for that purpose (rule 4). Additional 
work was needed to make diameter distribution predictions past 1940. The difference between 
the pre- and post-1940 periods was the number of diameter classes that could be compared (for 
pre-1940, all diameter classes 4 inches and greater could be compared while, after 1940, only 
diameter classes 8 inches and greater could be compared). 


The choice of what statistics and tests were appropriate for comparing predicted to the 
actual diameter distributions (rule 7) proved to be difficult. After reviewing or trying a 
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number of alternatives, none of them proved either applicable or useful (see appendix I for 
further details on tests examined or tried). A test not examined in this study but which 
could prove useful in future studies is described by Freese (1960). 


The lack of a useful statistical test, therefore, necessitated a subjective evaluation of 
the validation runs. This decision led to an additional problem: how to conveniently sum- 
marize, for reporting and comparative purposes, the large amount of detailed individual 
subplot validation information produced in validating for the uneven-aged plots. After 
consideration, I concluded that two statistics adequately reflected the average behavior 
exhibited by just the validation subplots. These statistics were computed by first averaging 
across all the plot's subplots the predicted and actual diameter distribution for each growth 
period, combining blackjack and yellow pine in the process. Next, the mean difference be- 
tween average actual and average predicted diameter distributions (the first statistic) and 
the variance of these differences (the second statistic) were computed across diameter 
classes for each plot and growth period. 


Use of averages incorporating both the model building and the validation data sets would 
appear to violate the fourth rule of validation (that is, the validation data should be indepen- 
dent of the model building data). But the apparent violation is not serious because much of 
the basic validation information produced and examined (but not reported here) for decision- 
making pertained to just the subplots from the validation data set. 


For the even-aged validation data, the mean difference across diameter classes, between 
the actual and predicted diameter distribution for each individual plot and growth period and 
the variance of these differences were used as summarization statistics. Statistics produced 
from average diameter distributions across all plots are meaningless for the even-aged data 
See. 


Results of Validation 
PREDICTED UPGROWTH; ACTUAL MORTALITY, CONVERSION, INGROWTH 


The first set of validation runs used predicted upgrowth and actual mortality, cutting, 
conversion, and ingrowth rates to assess the predictive capabilities of the various upgrowth 
models. 


In addition to the two blackjack pine basal area growth equations tested, two methods for 
applying the equations were also tested. One method used the average basal area growth for 
the diameter class to advance the trees, and the other divided the diameter class into three 
more classes and then used estimates for the lower, middle, and upper basal area growths of 
the diameter class to advance each third separately. The combination of basal area growth 
equations and advancement techniques resulted in the evaluation of four different upgrowth 
models. 


For simplicity, the actual mortality, conversion, and cutting rates were expressed as 
proportions rather than number of trees. This approach has the advantage of eliminating the 
problems associated with removals being larger than the number of trees in a diameter class. 
The method, however, exaggerates differences between predicted and actual upgrowth because 
any difference between the resulting predicted and actual diameter distributions will result 
in different numbers of trees being removed due to mortality, cutting, or conversion. 


Examining these runs revealed that, on the uneven-aged plots, the blackjack pine basal 
area growth equation developed with both even- and uneven-aged data (BJBAG1) behaved the same 
as the equation developed with just the uneven-aged data (BJBAG2). As a result, BJBAG1 was 
adopted as the final basal area growth equation for blackjack pine. 


In addition, the upgrowth rate appeared too high. This led to an immediate suspicion of 
the proposed correction for log bias because of the magnitude of that correction (a 35 percent 
increase for blackjack pine and 140 percent increase for yellow pine) and because the justifi- 
cation for that correction was not too firm. Therefore, I decided to try two additional 


correction schemes: no correction, and a correction of In(Y/Y), where Y is mean basal area 


growth and Y is mean predicted basal area growth with no correction for log bias. I chose 
the latter because the correction was about halfway between the original correction and zero 
correction. 


Six validation runs were made on each plot, each run representing a different log cor- 
rection and one of the two methods for applying the equations (that is, average basal area 
growth or basal area growth in thirds). From these runs I concluded that the best log bias 
correction was no correction, and the best method for applying the equations was to divide the 
diameter classes into thirds and project each third. 


The conclusion that a zero (when added to intercept terms of the log model) correction is 
best of the three methods should not be construed to mean that a zero correction is the best 
of all possible ones. A more likely interpretation is that the correction should be smaller 
than the two values tested. Obviously, the correction method posed in appendix D does not 
produce appropriate corrections for log bias. Until a more theoretical basis for correcting 
log bias in the nonnormal case is developed, the most reasonable approach may be to not correct 
the equations directly for log bias, but rather, to build in a calibrating routine that will 
correct for log bias, site differences, and so forth. 


Table 18 aids comparison of the effects of adding predictive components upon accuracy (as 
measured by average differences between predicted and actual number of trees in a diameter 
class) and precision (as measured by the variance of the differences). The first two columns 
of values provide the mean difference and the variance of differences for each plot and growth 
period (since initialization) for the set of runs using the final predicted upgrowth model and 
actual mortality, conversion, and ingrowth. If a model predicted upgrowth perfectly, the 
resulting values would be zero. While the values do differ from zero, their magnitudes indi- 
cate that upgrowth predictions are reasonably good. 


As expected, results also indicate that, as the number of growth periods from initial- 
ization increases, the mean differences and the variance of the differences also increase. 


PREDICTED UPGROWTH AND MORTALITY; ACTUAL CONVERSION AND INGROWTH 


I tested two blackjack pine mortality equations in this set of runs: one equation was 
developed using the uneven-aged data, and the other using both the even- and uneven-aged data. 
I developed the single yellow pine equation using uneven-aged data. 


Analysis of the runs showed the even- and uneven-aged blackjack pine equation evidently 
inferior to the uneven-aged blackjack pine equation for predicting uneven-aged mortality. 
The former equation also did not predict even-aged mortality well. This finding was not 
surprising because of the high level of snowbreak loss on the even-aged plots. 


A different approach for predicting even-aged mortality was tried. I reclassified the 
snowbreak mortality as a "cutting" loss (that is, not endemic), and predicted the remaining 
mortality using the uneven-aged mortality equations. The result was greatly improved runs 
over those using even- and uneven-aged mortality equations. 


This finding indicates that if catastrophic losses (such as snowbreak, fire, insects, or 
disease) are treated in the same fashion as cutting losses, endemic losses in even- and uneven- 
aged stands could be predicted using the same mortality equations, if the equations are 
developed appropriately. 


The third and fourth columns of table 18 provide the mean difference and the variance of 
differences for the runs using the final upgrowth and mortality models and actual conversion 
and ingrowth. A comparison with the values derived from the previously discussed set of runs 
indicates that adding the mortality functions has not greatly affected accuracy and precision. 
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Table 18.--Mean difference and variance of differences between average actual and average predicted number 
of trees per diameter class for a specified plot, growth period, and simulation type 


: Predicted upgrowth Predicted upgrowth, All components 
Plot Growth eee oe and mortality mortality, conversion predicted 
period Mean Variance of Mean Variance of Mean Variance of Mean Variance of 
difference difference difference difference difference difference difference difference 
61 1 0.0028 0.015i 0.0135 0.0153 0.0135 0.0153 0.0283 0.0138 
2 .0027 .0495 .0183 . 0486 .0183 . 0486 .0568 0371 
3 .0021 .0615 0231 .0640 .0231 .0640 .0683 .0563 
4 .0092 .0101 .0263 .0118 .0263 .0118 .0274 .0107 
5 .0082 .0102 0279 .0126 .0282 .0126 .0366 .0104 
6 .0124 .0244 .0268 40255 .0271 .0255 .0474 .0419 
¥ 5 .0126 .0327 -0295 .0366 .0300 .0365 20597 .0573 
8 .0261 .0395 .0450 .0459 .0455 .0457 .0800 .0713 
10 . 1000 - 4436 «1239 .4596 .1242 -4581 . 1600 . 4667 
62 1 .0000 .0093 .0080 .0094 .0080 .0094 -.0027 .0294 
2 -.0007 .0166 .0102 .0162 .0102 .0162 -.0050 .0511 
3 -.0010 . 0267 .0143 .0260 .0143 .0260 .0105 .0475 
4 .0334 .0101 .0492 .0092 .0492 .0092 .0516 .0106 
5 .0221 .0110 .0392 .0092 .0422 .0092 .0476 .0132 
6 .0250 .0214 .0455 .0178 .0461 .0178 20592 .0269 
if - 0366 .0541 0563 .0466 0571 .0465 .0700 .0541 
8 bly 35133 .1353 . 2966 lL Soi .2962 .1416 2529 
9 .0582 rly . 0808 .1610 .0813 . 1608 .0826 .1828 
71 1 .0003 .0292 - .0006 .0297 -.0006 .0297 . 0086 .0409 
2 -.0006 . 0806 .0003 .0844 .0003 .0844 -.0406 2727 
3 0627 -0335 .0645 .0323 .0645 .0323 .0727 .0423 
4 .0232 -0310 .0312 .0295 20312 .0295 .0497 .0260 
5 -.0343 .0426 -.0211 .0404 -.0243 .0407 -.0071 .0358 
6 -.0378 .0304 -.0228 0292 -.0244 .0295 -.0425 .0581 
8 -.2118 . 3379 -.1882 22975 -.1884 ZOU -.3484 1.6092 
72 1 . 0003 .0392 -.0017 .0363 -.0017 .0363 -.0042 .0588 
Zz .0005 .0710 .0038 .0692 .0038 .0692 -.0316 .2209 
3 . 0694 .0419 SO71S -0428 .0715 .0428 .0765 .0504 
4 .0466 .0403 .0534 .0422 0534 0422 .0631 .0481 
5 .0247 .0381 .0356 .0437 .0328 .0443 .0372 .0467 
6 .0811 .1905 AU / . 2096 .0938 =209 .0603 .0815 
8 -.0550 . 3290 -.0321 . 3246 -.0326 . 3256 -.2208 1.3616 
TW-5 1 .0000 39.87 .1744 40.67 -- -- -- -- 
2 -.0082 64.12 a2991 65.94 -- -- -.3419 40.49 
TW-6 i .0000 59.66 21567 60.86 -- -- -- =— 
2 .0100 78.00 . 1482 Si 271 -- -- -.5644 65.22 
TW-7 1 . 0000 20.96 .0622 21726 -- -- -- = 
2 .0000 20.87 0875 21.02 -- -- -.2291 9.80 
TW-10 i .0000 221.37 . 2786 226.01 -- -- == == 
2 .0000 164.43 -4778 ul7/esak -- -- -1.9250 51.86 
TW-13 if . 0000 39.24 e203 40.80 -- -- -- == 
2 . 0000 53.96 . 3400 56.99: -- -- -1.8774 32.58 
TW-17 if . 0000 20.74 .1022 21.14 -- -- =— == 
2 . 0000 Sieh bree 37.74 -- -- -.4680 35.55 


PREDICTED UPGROWTH, MORTALITY, AND CONVERSION; ACTUAL INGROWTH 


Conversion from blackjack pine to yellow pine is applied at the end of each fourth growth 
period. Therefore, only simulation runs longer than four growth periods will show the effect 
of adding predicted conversion upon the resulting predicted diameter distribution. For simu- 
lation runs (such as the even-aged runs) under five growth periods in duration, the results 
are identical to those runs using actual conversion and ingrowth. As a result, the even-aged 
data could not be used to check the predicted conversion models. 


The fifth and sixth columns of table 18 present the mean and variance of differences for 
the simulator using the final predicted upgrowth, mortality, and conversion models, and actual 
ingrowth. Results indicate that using predicted conversion does not reduce accuracy and 
precision by much. 


ALL COMPONENTS PREDICTED 


The final component, ingrowth, was developed using the uneven-aged data alone and, for 
this reason, is not applicable to the even-aged validation plots. The addition of the ingrowth 
component reduces both accuracy and precision (see the last two columns of table 18), which 
is particularly noticeable on the longer runs and on the managed, uneven-aged plots (plots 
7A cami y/:2i)i. 


Because an objective was to produce a simulator useful in predicting even-aged stands as 
well, a substitute for the ingrowth model was devised and then tested on the even-aged valida- 
tion plots. This consisted of extending the use of the upgrowth and mortality models into the 
1-, 2-, and 3-inch diameter classes as well. If the number of trees existing below the 1l-inch 
diameter class and the possibility of additional regeneration are both insignificant for the 
stand, then ingrowth should consist of those trees from the 1-, 2-, and 3-inch classes surviving 
and upgrowing into the 4-inch or larger diameter classes. Therefore, if the available upgrowth 
and mortality equations can be extended into these classes, reasonable ingrowth predictions 
should result. 


The values listed in the last two columns of table 18 for the Taylor Woods plots are the 
results of applying the proposed technique for two growth periods. On all plots, the technique 
both increased the mean difference per diameter class and reduced the variance of the differences. 
The resulting mean differences, however, are still small when compared to values for average 
actual number of trees per diameter class that ranged from 9.38 (for TW-7) to 77.90 (for TW- 

10). I concluded from these findings that the technique is an adequate substitute for an 
ingrowth model if it is applied to a stand with an average diameter of at least 4 inches. 
This is necessary because all of the Taylor Woods plots fell in this category. 


Predicting Subplots and Averaging vs. Predicting Plot Averages 


Of those uneven-aged whole stand simulators developed and published, all have used 
relatively small plots as their data base. Moser (1972) used 1/5-acre plots. Later Moser 
(1974) used plots ranging in size from 1 to 4 acres. Finally, Ek (1974) used 1/7-acre plots. 
The apparent assumption is that these small plots have the same structure and dynamics as the 
total stand. This assumption is probably met in truly all-aged stands; however, if clumping 
exists (as it does in ponderosa pine), the assumption probably cannot be met. 


To test this assumption, the average plot diameter distributions were predicted and 
compared to the values previously obtained by averaging the predicted subplots. Table 19 
presents mean differences and the variance of the differences for each plot and growth period 
using the method of predicting plot averages. Comparing these values to those found in the 
last two columns of table 18, I concluded that predicting plot averages generally reduces 
both accuracy and precision. The loss, however, appears not too severe for most potential 
uses of a whole stand simulator. The conclusion is fortunate since using the simulator for 
answering managerial questions would be greatly restricted if it were necessary to first 
simulate subplots and then average to get stand values. 
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Table 19.--Mean difference and variance of differences between actual and predicted average 
plot number of trees per diameter class for a specified plot and growth period 


Plot Growth Mean Variance of 
period difference difference 
61 1 0.0275 0.0226 
2 .0563 .0595 
3 . 0636 .0939 
4 .0474 20025 
5 0742 .0168 
6 . 1076 .0746 
7 . 1368 .1100 
8 GAS eA Sal 
10 2A ESQ IT 
62 1 =. 0224 .0449 
2. -.0464 .0647 
5 -.0607 SS 4 
4 One .0070 
5 .0263 .0057 
6 . 0284 .0134 
7 . 0300 .0119 
8 .0642 . 0664 
9 -.0453 P2520 
71 1 -.0137 .0780 
2 -.0847 .4184 
3 .0564 .0426 
4 0362 .0283 
5 .0189 .0292 
6 -.0486 .0807 
8 -.4189 2.0891 
72 i -.0264 .0991 
2 -.0773 3569 
3 .0615 .0399 
4 0409 .0368 
5 . 0469 .0196 
6 .0378 50223 
8 -.3097 Le7S29 


RESULTS OF STUDY 


Analysis of the validation runs demonstrates that the simulator provides reasonable results 
for runs at least 40 to 50 years in duration (that is, 8 to 10 growth periods). This is 
fortunate because cutting cycles can be from 5 to 30 years long, and the method for determining 
optimal diameter distribution requires good estimates of the future stand at the end of the 
cutting cycle. In addition, solving the optimal conversion strategy problem could require 
simulation runs as long as two or three cutting cycle lengths. 


The final equations that make up the simulator are found in table 20. Those equations 
incorporating predicted basal area growth as an independent variable (the mortality and 
through-growth equations) have had their regression coefficients changed to adjust for the 
final removal of the log bias correction. A description of the control cards for the simulator 
can be found in appendix J. 
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Table 20.--Summary of final equations for simulating uneven-aged stand dynamics 
in ponderosa pine/Arizona fescue habitat type 


BASAL AREA GROWTH EQUATIONS 


Average Basal Area Growth 
Blackjack pine 
In(BJBAG) = -8.51836897 + 1.16754330(in(D)) - 4.00970143E-02(D) 
-3.84298771E-03(LBA») - 7.15483662E-03(MBA>) - 1.58234269E-02(UBA,) 
-3.26097273E-01(A,*1n(D)) + 8.80676713E-01(A3) + 1.0(1in(S)) 


BJBAG = EXP(in(BJBAG) ) 
Yellow pine 
In(YPBAG) = -15.2464932 + 4.27656656(i1n(D)) - 1.8316162E-01(D) 
-7.27361567E-05(LBA>)* - 9.11165626E-04 (MBA .) 2 
-2.41462106E-04(UBA,) 2 - 1.05776062(A,) + 1.0(in(S)) 


YPBAG = EXP(in(YPBAG) ) 


Fast, Moderate, and Slow Basal Area Growths 
Blackjack pine 
Fast BJBAG = EXP(0.56557130 + 0.12995917(A>) + In(BJBAG) ) 
Moderate BJBAG = EXP(0.09120012 + 0.12946423(A,) + In(BJBAG)) 
Slow BJBAG = EXP(-0.65677142 - 0.25942340(A,) + In(BJBAG)) 


Yellow pine 
Fast YPBAG = EXP(0.70338067 + 0.7280557(A3) + 1n(YPBAG)) 


Moderate YPBAG = EXP(0.24266226 + 0.59168455(A3) + In(YPBAG)) 
Slow YPBAG = EXP(-0.94604293 - 1.31974025(A3) + In(YPBAG)) 


MORTALITY EQUATIONS -- Proportion of trees dying in next 5-year period 
col 
PM = (1.0 + EXP(-X)) 


Blackjack pine 
X = -5.372649 + 1.270862E-03(D2) - 12.96478 (A, *BJBAG) 


Yellow pine 
X = -3.612474 - 4.226903(YPBAG) - 22.24186(A,*YPBAG) + 5.597303E-04 (D2) 


CONVERSION EQUATION -- Proportion of blackjack pine converting to yellow pin 
in each fourth growth period ; 


il 
PG = (1.20: +, EXP'\(-¥)) 
Y = -3.777269 + 0.5012605(D) - 6.59746SE-03(D2) - 0.0982176(S) 


UNEVEN-AGED INGROWTH EQUATIONS 
Total Ingrowth 


TI = boS(BACL,)?1(1.0 + BACL>)?2EXP(b3BA?* + be (A,*BA?4)) 


bo = 0.48630258 b3 = -3.97721395E-06 
b, = 0.6594918 by = 2.7667388 
bop = -0.49280582 bs = -1.33299951E-06 
Through-Growth -- Proportion of total ingrowth growing into 5-inch diameter class 


3 
PT = (1.0 + EXP(-Z)) 
Z = -2.636336 + 38.80808(BJBAG-4") - 0.03419626(A,*TI) 


(continued) 
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Table 20 (continued) 
Ingrowth into 4-inch diameter class 
ING ,,, =a aie) (Glia) — te) 
Ingrowth into 5-inch diameter class 


ING.,, = UAT 


HEIGHT EQUATIONS 
Even-aged 
aera a z Daves 
Hp = MH - 0.88015037(MH - 4.5)(1.0 - MD. 


Uneven-aged 


-2 
Hy = G9) + ey so2 . EUS) ; 


Co = 4.5 

ce) = 13.178649 
co = 0.71631005 
c3 = -4221.6528 


DEFINITION OF VARIABLES 
BJBAG = Predicted average 5-year basal area growth, in square feet, of blackjack 
pine 
YPBAG = Predicted average 5-year basal area growth, in square feet, of yellow 
pine 
BJBAG-4"' = Predicted average 5-year basal area growth in the 4-inch diameter 
class of blackjack pine 


PM = Proportion of trees dying in the next 5-year period 


PC = Proportion of blackjack pine converting to yellow pine in each 
fourth growth period 


TI = Total ingrowth for uneven-aged stands 


PT = Proportion of total ingrowth growing into the 5-inch diameter 
class for uneven-aged stands 


ING,,, = Ingrowth into the 4-inch diameter class for uneven-aged stands 
ING.,, = Ingrowth into the 5-inch diameter class for uneven-aged stands 
Hh = Average height of the Dth diameter class 
D = Diameter class size 
S = Minor's (1964) site index 


MBA» = Basal area in the given diameter class plus the two adjoining 
larger and smaller diameter classes 


LBAz = Total basal area below the smallest diameter class in MBA» 
UBAy = Total basal area above the largest diameter class in MBA> 
BA = Total basal area 
= MBA, + LBA» + UBAy, for a specified diameter class 


TIME = Number of 5-year growth periods since last cutting 
DM = Maximum diameter class size that an even-aged stand has achieved in its development 
MH = Maximum height as predicted by Minor's (1964) equation 
S - 1.4003(VA - 10) + 0.1559(S) (VA - 10) 
A = Breast height stand age, 20<A<140 


Basal area in the 4- through 6-inch diameter classes 
Total basal area above the 6-inch diameter class 
-0.244178 + 1.244178*EXP(-(1.176471 - 0.019607*TIME) 3) 
-0.095336 + 1.095336*EXP(-(1.25 - 0.020833*TIME) *) 
-0.000203 + 1.000203*EXP(-(1.428571 - 0.023809*TIME) ©) 
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CONCLUSIONS 


A simulator for predicting uneven-aged stand development was constructed and validated 
for the ponderosa pine/Arizona fescue habitat type of the Southwest. The structure and dynamics 
of the simulator make it potentially useful for answering even- and uneven-aged management 
questions. Some results and conclusions formulated during the construction and validation 
process are: 


1. The distribution of number of trees within l-inch diameter classes is uniform. 


2. Diameter class basal area growth rates of even- and uneven-aged stands are logically 
interrelated if the structure of within stand competition is incorporated in the model. 


3. No reasonable log bias correction could be found for the log of basal area growth 
equations when residuals are not normally distributed. In the case of nonnormality, the 
development of a calibrating routine is recommended as an alternative to a log bias correction. 


4. Even- and uneven-aged endemic mortality rates can be logically interrelated if cata- 
strophic mortality is treated in a fashion similar to cutting mortality. 


5. The ingrowth component is the weakest link in the simulator. The development of a 
regeneration, seedling, and sapling component is highly desirable. 


6. The use of the chi-square statistic for testing "goodness-of-fit" in validation is 
highly questionable. 


7. The use of regression statistics for testing "goodness-of-fit" in validation did not 
prove to be very valuable. 


8. The simulation of 2.5-acre subplot values followed by the calculation of their average 
to obtain mean plot values produced different results than simulating plot averages directly. 
Because all published uneven-aged stand simulators have used small plots, this result could 
prove a problem for those interested in simulating plot averages. For the simulator developed 
in this study, however, analysis also revealed that the difference between the two methods was 
not too severe. 
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APPENDIX A 
DATA PROBLEMS 


The even- and uneven-aged plots used in this study provide an unusual modeling opportunity 
for western species. Unfortunately, their use also posed numerous modeling problems: 


1. Narrow geographic and site range.--All of the uneven-aged plots are situated within a 
square mile tract, and the even-aged plots are about 1 mile away. The uneven-aged plots have 


a narrow range in site and, because of this, site index as an independent variable may not be 
significant for predictive purposes. In addition, the even-aged plots have an average site 
index higher than uneven-aged plots. While the differences in site index may be insignificant 
for the uneven-aged plots, the difference between the uneven- and even-aged plots may well be 
Significant. All this could have contributed to the problems encountered while trying to 
incorporate site index in the various component models of the simulator. 


Fortunately, one mitigating factor for these problems is that all the plots (both even- 
and uneven-aged) occur on the same habitat type, Pinus ponderosa/Festuca arizonica (ponderosa 
pine/Arizona fescue), and it is believed that the plots and their site indices are representa- 
tive of the type. Therefore, the simulator developed in this study is probably applicable to 
that habitat type, at least in northern Arizona. 


2. Changes in diameter measurement limits.--Ideally, the smaller the lower diameter 
limit, the better. Large, lower diameter limits ignore too much of the stand that is important 
in understanding and modeling stand dynamics. The changing of lower limits on the Fort Valley 
uneven-aged plots caused problems and raised several questions. A lower diameter limit of 7.6 
inches would be necessary to utilize the full-time range of the data. This would provide six 
to nine 5-year growth periods of data and would incorporate cuttings at several times. It 
would, however, also ignore about one important third of the stand--all trees under 7.6 
inches. 


Conversely, a lower diameter limit of 3.6 inches would restrict the data to only the 
early time periods. This would provide three to four 5-year growth periods, but it would also 
eliminate many of the cuttings. 1 decided to use only the early data because: (a) a larger 
portion of the stand would be modeled, and (b) by including more of the total stand, the 
possibility of meeting objective number 2 would be increased. 


Because of the lower diameter limit of 3.6 inches for the uneven-aged stands, all data on 
the even-aged plots below 3.6 inches were ignored. Several plots in the Taylor Woods study 
have many small trees and, as a result, much of their structure was eliminated. A greater 
portion of the number of trees in the total stand is usually located ani the 020-tors.5-imeh 
class of uneven-aged stands than in most even-aged stands, with the exception of very young 
even-aged stands. While a lower diameter limit of 3.6 inches is an improvement over 7.6 
inches, it is still far from the ideal solution of modeling the total stand. 


3. Intrusion upon uneven-aged plots.--Unfortunately, a major paved highway runs through 
pieces of all the uneven-aged plots, and this could affect their stand dynamics. Not only 


would these disturbed subplots have a serious "edge effect,'' but the two remaining subplot 
pieces would not have the same competitive interaction as occurs on the undisturbed subplots. 
The number of subplots seriously affected by the highway are given in table 21. 
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Table 21.--Number of subplots by size, plot number, and presence of highway 


Subplot Plot 6l Plot 62 Plot 71 Plot 72 

size No No No No 

in highway Highway highway Highway highway Highway highway Highway 
acres intrusion intrusion intrusion intrusion intrusion intrusion intrusion intrusion 


0.50 3 
.70 if 


a oo 
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.50 20 28 1 23 


AWW MM MDM MN DN NN NH KH 
Co ‘= 
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25 9 Sill 1 24 8 Sill 1 


Number 

of 

subplots 

seriously 

affected 

by 

highway 6 0 5 0 


4. Variable plot sizes.--As mentioned earlier, sizes on the uneven-aged subplots range 
from 0.5 to 3.1 acres (table 21), and the even-aged plots from 0.75 to 1.25 acres. Curtis and 
Pope* hypothesize that '"'small" plots can result in erratic estimates of stand attributes, such 
as basal area, because of within-stand clumping. As a result, the deviation of plot estimates 
of growth and basal areas from their respective stand means will be correlated and, therefore, 
the regression of diameter class growth as a function of basal area will also reflect plot 
size. They conclude that this problem would not produce biased results if the regression 
equation were applied to inventory plots of the same size and design, and if the stands to 
which they are applied had similar spatial patterns. Recognizing that these conditions are 
seldom met, Curtis and Pope suggested that the problems would be minimized if large plots were 
used in developing growth equations. 


Plot size is important to this study because the two data sets (even- and uneven-aged) 
have different plot sizes, both within and between the data sets. based on the Curtis and 
Pope argument, the coefficients for the "separate" models (even- and uneven-aged) may be con- 
founded by the plot size differences, making a statistical combination more difficult. To 
help minimize the problem with the uneven-aged subplots, all subplots below 2 acres were 
eliminated. This action, coupled with the realization that the even-aged plots are fairly 
homogeneous (because of stand treatment) and good sized, might reduce the problem of within- 
stand clustering and the possibility of locating a subplot within one cluster of the stand. 


2Curtis, Robert O., and Robert B. Pope. 1972. Some considerations in design of growth 
studies and associated inventories of the western National Forests. Draft of the internal 
report. USDA For. Serv., Pac. Northwest For. and Range Exp. Stn., Portland, Oreg. 


47 


5. No height-growth information on uneven-aged plots.--The lack of height-growth data on 


the uneven-aged plots did preclude the development of a common even- and uneven-aged height- 
growth model, but its absence was not critical to the dynamics of the basic components in the 
simulator. Because some uneven-aged height-diameter data were available from both the School 
of Forestry at Northern Arizona University and the Rocky Mountain Forest and Range Experiment 
Station at Flagstaff, I decided to develop separate height-diameter relationships for even- 
and uneven-aged stands. While height-diameter relationships are not ideal substitutes for a 
height-growth model because they do not incorporate the effect of stand density upon height, 
their use does provide a means of assessing product potential. 


The School of Forestry data came from student class measurements of l-acre square subplots 


made on plot 61 during 1964 and 1967. The Rocky Mountain Forest and Range Experiment Station 
data came from one-time measurements on selected subplots of plots 71 and 72. 
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APPENDIX B 
DISTRIBUTION OF TREES WITHIN A DIAMETER CLASS 


The uniformity of tree distribution within a diameter class is likely to be influenced by 
the combination of diameter class size and the number of trees in the diameter class. The 
fewer the number of trees in a class, the more difficult it would be to prove nonuniformity; 
and the number of trees in a class is, in part, determined by the diameter class size. 
Therefore, the smaller the diameter class size, the higher the probability that the assumption 
of uniformity can be met. 


To test for uniformity, about one-fourth of the subplots reserved for model building were 
randomly selected ffom both the uneven- and even-aged plots. All applicable measurement 
periods on the selected subplots were examined, and those l-inch diameter classes with five or 
more trees in them were then chosen for testing. 


A chi-square "goodness-of-fit" test was selected because it is appropriate for testing 
the distribution of a discrete random variable (the number of trees in a class). For diameter 
classes with five to nine trees in them, the diameter class was divided into five equal sub- 
classes, and for diameter classes with 10 or more trees, the number of subclasses was 10. 

The expected number of trees in each subclass was computed by dividing the total number of 
trees in the class by the appropriate number of subclasses. 


The decision to use two sizes of subclasses was based on the suggestion by Snedecor and 
Cochran (1967) that the expected value in any subclass should be greater than or equal to one. 
Since the number of trees in a diameter class is usually low in the Southwest, the inclusion 
of the five-to-nine size group spreads the testing over more diameter class sizes. For example, 
of the 2,195 diameter classes with trees examined in this phase, 1,481 had less than five 
trees, 417 had five to nine trees, and the remaining 297 diameter classes had 10 trees or more 
(table 22). 


Table 22.--Results of testing within diameter class diggribution for each plot by number of 
trees in diameter class e 


Number of diameter classes 


Plot 5 to 9 trees 10+ trees 
1 to 4 trees Uniform Nonuniform Uniform Nonuniform 

61 499 147 1 39 2 

62 415 79 0 8 0) 

71 257 98 0 94 2 

U2 299 87 0) 101 0 
Taylor Woods ll 5 0 49 2 

Total 1,481 416 1 291 6 


Lestimpuwasi done at the 99 percent) level of significance im order not to reject the null 
hypothesis of a uniform distribution unless marked deviations were found. The degrees of 
freedom for the test were computed as the number of subclasses minus one because it was not 
necessary to estimate distribution parameters. Results of the tests are found in table 22. 


Of the 714 diameter classes tested, only seven were not uniformly distributed. While it 
was earlier hypothesized that the within diameter class distribution might depend upon number 
of trees in the class and perhaps stand structure, the results indicate that these factors are 
not as important as first thought. For example, on the even-aged plots tested, 20 of the 21 
diameter classes with 50 or more trees in them were uniformly distributed, indicating that the 
within diameter class distribution is not closely correlated with the number of trees or 
structure of the stand. 
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APPENDIX C 


DEVELOPMENT OF DIAMETER GROWTH MODELS 
Prior Findings 


Diameter growth is influenced by several factors including productivity. A measure of 
productivity is site index (the height obtainable by free growing trees at a base age). 
Minor's (1964) site index was developed for even-aged stands in northern Arizona, and its use 
is recommended for the area.*? The application of site index to uneven-aged stands is compli- 
cated by the difficulty of finding suitable free growing site trees. This complication is 
reduced to some extent in the uneven-aged stands because they have relatively low stocking 
levels; however, this problem still causes difficulties. Another problem with site in the 
study area is the previously mentioned narrow range of site indices computed for the plots, 
which may make the inclusion of a meaningful site measure difficult. 


Another factor related to productivity that also influences diameter growth is rainfall. 
Fritts and others (1965) found that diameter growth in southwestern ponderosa pine is influenced 
by precipitation and temperature of the autumn, winter, and spring prior to the growth period. 
Pearson and Wadsworth (1941) found that gross 5-year volume increment appeared to be related 
to the amount of precipitation in the period. Finally, Pearson (1936) hypothesized that a 
fall off in 5-year diameter growth in the second period after thinning was due to lower precipi- 
tation in that period. 


Total stand density can also influence diameter growth. For uneven-aged stands in the 
Southwest, Krauch (1940) and Pearson (1950) both reported that a reduction in stocking through 
cutting increased diameter growth. After evaluating growth in young, even-aged clumps of 
uneven-aged stands, Cooper (1961) found that average tree diameter was "'...largely determined 
by stand density."’ For even-aged ponderosa pine stands in the Southwest, the effect of density 
upon growth has been widely reported (Pearson 1936; Hornibrook 1936; Krauch 1949b; Gaines and 
Kotok 1954; Cooper 1960; Larson and Minor 1968; Schubert 1971, 1974; and Myers and others 
1976). In all cases, an increase in density reduced growth. 


The correlation of diameter class size with growth has also been reported in the Southwest. 
Some have found that, following cutting in uneven-aged stands, diameter growth was not as 
closely correlated with diameter class size as might be expected (Pearson 1933, 1950; Pearson 
and Marsh 1935). Others, however, have found that diameter growth did vary markedly between 
classes in cut stands (Krauch 1937) and in uncut stands (Pearson 1950). For even-aged stands, 
both Pearson (1936) and Krauch (1949a, 1949b) reported differences in average diameter growth 
between diameter classes. 


One feature mentioned frequently in the literature is the difference in growth rates 
between "'yellow'' pine and ''blackjack"' pine. In all cases, diameter growth was greater in 
blackjack pine than in yellow pine (Pearson 1933, 1950; Pearson and Marsh 1935; and Krauch 
1934),: 1937)\- 


While not mentioned specifically in the literature for the Southwest, one additional 
factor that does affect diameter growth is the position of the diameter class in the stand 
(that is, how much competition a given stand receives from above and below it). A measure of 
this type can be thought of as an index of the competitive structure, and it is intimately 
related to stand structure. A competitive structure index might, therefore, explain differences 
between even-aged and uneven-aged growth rates. Competitive structure was one of the critical 
factors for meeting objective number two of this study. 


3Personal communication with Gilbert H. Schubert, principal silvilculturist with USDA For. 
Serv., Rocky Mountain Forest and Range Experiment Station, Flagstaff, Ariz.; retired. 


General Model Form 


When modeling diameter growth, one of the first decisions is the general form that the 
model should take. The three most common forms suggested for growth modeling are the following: 


Y=b'iex+e (1) 

Yes flax) Sc (2) 

Y= 20;<)\- © (3) 
where 


Y = the dependent growth variable 


an array of independent variables that may themselves be functions of basic growth 
factors 


Ids 
tI 


a, b = arrays of model parameters to be estimated in the appropriate fashions 
f(,) = a nonlinear function 
€ = a random error. 


Model form (1) assumes that growth is a linear, additive function of transformed independent 
variables. The estimation of the parameters of the model is done through ordinary least 
Squares regression techniques. Vuokila (1965) used this form when he modeled percent diameter 
growth of individual trees. 


Model form (2) assumes that growth is a nonlinear function of the independent variables. 
Because the error term, €, is additive, the appropriate technique for estimating the model 
parameters is nonlinear, least squares regression (Kmenta 1971). 


Model form (3) differs from (2) in how the error term is introduced. In model (3), the 
error term is multiplicative instead of additive. This allows the error structure model to be 
linearized through the use of natural logarithms. For this model form to work, it is necessary 
that the parameters (b) of the nonlinear function can also be linearized through the log 
transforming process. The resulting model is of the form: 


Lal AY Wo D8 0) G:(a0)' or 'eX (4) 
where 


g(x) = an array of appropriately transformed independent variables to linearize the 
parameters of the the function f(b,x) 


eae aliri((e) 


If these assumptions are met, then estimation of the parameters in model (4) can be done 
by ordinary least squares regression techniques. Model form (4) has been used by Lemmon and 
Schumacher (1962) to model individual tree periodic radial increment, and by Cole and Stage 
(1972) and Stage (1973) to model individual tree periodic basal area increment. 


Another choice faced by the modeler of diameter growth is the form of the dependent 
variable. Following are some of the choices: radial growth, diameter growth, basal area 
growth, or any of the previous three expressed as a percentage. Cole and Stage (1972) and 
Stage (1973) settled upon the usage of basal area growth for two reasons: 
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1. Basal area growth is often nearly linear over short time periods, and this makes 
extrapolation of growth rates easier for growth periods different from that originally used in 
equation development. 


2. The log of basal area growth is often found to be more normally distributed with 
homogeneous variance than other dependent variables. 


The latter finding may also indicate that the residuals of the log of basal area growth 
are additive. If the assumption of normality, homogeneity of variance and additive, indepen- 
dent residuals can be met, then the ordinary least squares estimators, b, are the maximum 
likelihood estimators and the UMVUE's (Uniformly Minimum-Variance Unbiased Estimator), and 
various procedures for testing the significance of the model and its parameters can be validly 
applied (Kmenta 1971; Draper and Smith 1966). 


Given these various choices of model forms and dependent variables, I decided that the 
log of basal area growth would be used because: 

1. Previous experiences with modeling diameter growth in ponderosa pine* indicate that a 
nonlinear, multiplicative model best represents the interaction of the independent variables 
with themselves and their effect upon diameter growth. 


2. Cole and Stage (1972) and Stage (1973) have shown that basal area growth is often 
lognormally distributed with a multiplicative error term. This fits the requirements of model 


Gir 
The Random Error Component 


The two approaches for handling the random error component are stochastic modeling and 
deterministic modeling. In stochastic modeling, each component of the model with an error 
element is randomly assigned an error value from the proper distribution. The prediction from 
this component is then modified by this random disturbance, and this "randomized" prediction 
interacts with the other model components to produce a randomized estimate from the model 
(usually with an unknown error structure). This process is then repeated a number of times, 
each time using a new set of random errors, and the results of these numerous trials can be 
averaged to get an "expected" model estimate. This repetitive stochastic process is often 
called a "Monte Carlo" process. One disadvantage of this approach is the sometimes large 
number of computations needed to find the "expected" estimates. 


The deterministic approach assigns each component of the model its expected value. There 
is no random element. Each component interacts with other components in the appropriate 
fashion to produce model estimates. The advantage of deterministic models is that "expected" 
estimates are provided directly. 


Stage (1973) describes and uses a method in his prognosis model that maintains the relative 
Simplicity of the deterministic approach while introducing a stochastic element. He calls the 
approach a Monte Carlo "'swindle."" He says, ''The purpose is to produce a prognosis that over- 
all is the result of averaging many replications of the random process without actually having 
to carry out the replication.'' The random element of this "swindle'' is applied to the log of 
basal area growth equations in one of two ways. The particular method is dependent upon the 
number of trees sampled in the stand. 


“David W. Hann. Diameter growth equations developed for Mora County, New Mexico, and the 
Black Hills, South Dakota, inventories, on file with the Forest and Range Resources Evaluation 
Project, USDA Forest Service, Intermountain Forest and Range Experiment Station, Ogden, Utah. 
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For large samples, Stage randomly picks a deviate from the normal distribution of the 
residuals about the log of basal area growth regression equation and adds it to the estimate 
of log of basal area growth. The underlying assumption is that, with a large sample, the 
effects of the individual random deviations will average out over the stand. 


For small samples, each sample tree record is divided into three sample tree records. 
The number of trees in each record is a fixed proportion of the original number of trees 
represented by the old sample tree record. The proportion breakdown Stage used was 15, 60, 
and 25 percent, based on previous findings that an average stand had 15 percent suppressed 
trees and 25 percent dominant trees.° Each of the new sample tree records is assigned a 
growth rate computed by taking the average growth rate (as predicted by the log basal area 
growth equations) and adding a random component. For the 25 percent dominant trees, the 
random component is the expected residual value of the largest 25 percent of the normally 
distributed residuals about regression; for the 15 percent suppressed trees, it is the expected 
residual value of the lowest 15 percent of the normally distributed residuals, and similarly 
for the middle trees. As a result, the weighted log of basal area growth still sums to the 
average log of basal area growth, as predicted by the equation. At each simulation period, 
the sample tree records are split again until enough sample trees exist so that the first 
method can be used. 


Certain aspects of the second method seem applicable to this study. Within a diameter 
class, individual tree growth could be quite variable. By using only average diameter growth, 
all trees in the class will be assigned the same growth rate and, therefore, advancement to 
larger diameter classes will be identical. If, however, the number of trees in a diameter 
class were divided into thirds (each third representing the fast, slow, and moderate growers), 
then a separate growth rate could be assigned each third in the same fashion as used by Stage. 
If the three growth rates differed enough, this would then allow the trees in a diameter class 
to move into a wider range of larger diameter classes. If the distribution of residuals 
remains constant over time, this approach should provide a more realistic estimate of the 
number of trees changing diameter classes. Thirds were used in this study instead of Stage's 
proportions because the proportion of suppressed and dominant trees for southwestern ponderosa 
pine is unknown. 


Definition of Independent and Dependent Variables 


The previously discussed factors influencing diameter growth include diameter class size, 
productivity, rainfall, total stand density, vigor, and position of the diameter class within 
the stand. 


Over diameter class size (D), it is expected that basal area growth will climb to some 
peak value and then drop off. Where and at what magnitude the peak will occur is unknown. 
One function that would allow a wide range of peaking forms is the Weibull function. The 
Weibull function can be expressed as (Bailey and Dell 1973): 


c 
| 


e-] ec 
. & (5) * EXP (5) ,» where b and c are parameters 


or 


“Based on a presentation given by Albert R. Stage, principal mensurationist with the 
Intermountain Forest and Range Experiment Station, USDA Forest Service, in a spring 1976 
graduate seminar at the University of Washington. 
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Linearizing this function gives: 
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Generalizing this function further to provide for an even wider range of potential model forms: 


Cc bi 50 
In(y)ossbou Dylniexetaboxe, Gb>= ows 


Therefore, the inclusion of In(D) and D° in the log of basal area growth model would 
allow for a wider range of "peaking" functions. The values of cxuchosen werewiy WeSe]25, 225s 
3, 3.5, and 4. These values cover a wide range of forms from the Weibull function (Bailey and 
Dell 1973). 


The measure of productivity in this study is Minor's (1964) site index. It is expected 
that, when site index (S) is zero, growth would be zero, and, as site increases, growth would 
also increase. To model this effect, the proposed independent variable is the log of site 
index, in other words, In(S). This choice of independent variable is supported by the findings 
of Cole and: Stage (1972). 


Because of the findings of Fritts and others (1965) and Pearson and Wadsworth (1941), I 
decided to express precipitation in two fashions: average annual 5-year growth period rainfall 
(GRF), and average annual 5-year rainfall (ARF). Growth period rainfall is defined as the 
rain that fell in the interval from September through May previous to the given growing period. 
Based on a reasoning similar to that for site index, I decided that the log of GRF and the log 
of ARF were the appropriate independent variables. 


There are several different expressions of total stand density, including the following: 
total stand basal area per acre (BA), total number of trees per acre (T), and crown competition 
factor (CCF) (Larson and Minor 1968). In all cases, growth should maximize when the stand 
density values are near zero, and growth should decline as the values increase. An appropriate 
form for modeling this behavior is: 


Ya ey 
where 
y = predicted growth 
x = measure of stand density 


bg, b, = model parameters 


d = a power on x 


Linearizing this form produces the following independent variables: BAG, Tas and Cena 


Three potentially suitable values of d are 2, 3, and 4 (d could take on any positive value). 
These particular values were chosen because any one of them would provide for a slow decline 
in y when x is near zero, which is also an expected effect (that is, low stand densities 
will not influence growth). 


The vigor differences between blackjack pine and yellow pine were handled by simply 
modeling each of them separately. 
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The final class of independent variables is those indicating the position of the diameter 
class within the stand. For this purpose, Stage (1973) developed the independent variable 
"percentile in the basal area distribution" (PCT). Stage defined this variable as the basal 
area in all trees equal to or smaller in diameter than a given tree, divided by total stand 
basal area. In this study, PCT was defined in a similar fashion, substituting ''diameter 
class" for "tree.'’ Transforms of this variable include the following: I1n(PCT), PCT (Stage 


AZ 
1973), and PCT”. 


Another "position" variable that Cole and Stage (1972) and Stage (1973) have used is 
diameter of the tree divided by the diameter of the tree of mean basal area (DCP). This 
independent variable can also adapt to the diameter class situation by dividing diameter class 
size by the diameter of the tree of mean basal area. Potentially, useful transforms of this 
variable include DCP and in(DCP). 


The final method of representing the position of the diameter class within the stand 
consists of separating total stand basal area into three components. The middle basal area 
component is defined as that falling in a specified diameter class (MBA); or within the 
specified diameter class plus the adjoining larger and smaller diameter classes, if they 
exist (MBA>); or within the specified diameter class plus the adjoining two larger and two 
smaller diameter classes, if they exist (MBA3). The lower basal area component (LBA), LBA2, 
or LBA3) is the total basal area existing below the smallest diameter class in MBA |, and the 


upper basal area component (UBA;, UBAj, or UBA3) is the total basal area existing above the 
largest diameter class in MBA. . 


In this fashion, three sets containing three independent variables were derived: LBA, 
MBA,, and UBA,|; LBA», MBA>, and UBAz; and LBA3, MBA3, and UBA3. These variables should behave 
in the same fashion as the total stand density variables, and therefore the only additional 
transformation used was to also square each value. This resulted in six sets, each containing 
three independent variables. For convenience, the definition of all independent variables and 
their abbreviations can be found in table 23. 


Due to the desire to divide the residuals about the final models into fast, average, and 
slow growers, the appropriate dependent variable is the log of basal area growth for each 
tree. While the dependent variable is an individual tree value, it must be remembered that 
all independent variables are based only on diameter class information. 


Development of Equations 


Equations were developed for both blackjack pine and yellow pine using the least square 
regression program REX (Grosenbaugh 1967). REX is a powerful screening tool because it is 
designed to fit all combinations of independent variables following specified combinatorial 
rules. The screening statistic used is the relative mean square residual (RMSQR), which is 
defined as the mean square error about regression divided by the variance of the dependent 
variable. Therefore, a perfect fit results in an RMSQR value of zero, while an RMSQR value of 
one would indicate that the regression equation is no better at reducing squared residuals 
about regression than a simple mean. The advantage of the RMSQR over the more commonly used 
coefficient of determination (R*) is that the former reflects the reduction in the degrees of 
freedom caused by adding another independent variable. Therefore, it is easier to compare 
models with a different number of independent variables. 
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Table 23.--Definition of independent variables and abbreviations 


Abbreviation Independent variable 
D Diameter class size (e.g., for the diameter class 8.6 to 9.5, D = 9) 
S Site index 
GRF Average annual 5-year growth period rainfall 
ARF Average annual 5-year rainfall 
BA Total stand basal area per acre 
CCF Crown competition factor 
ie Total number of trees per acre 
PGT Percentile in the basal area distribution 
DCP Diameter class size divided by the diameter of mean basal area 
LBA} Total basal area under a given diameter class 
MBA} Total basal area in a given diameter class 
UBA, Total basal area over a given diameter class 
LBA» Total basal area under (a given diameter class - 1) 
MBA» Total basal area in (a given diameter class + 1) 
UBA» Total basal area over (a given diameter class + 1) 
LBA3 Total basal area under (a given diameter class - 2) 
MBA3 Total basal area in (a given diameter class + 2) 
UBA3 Total basal area over (a given diameter class + 2) 
TIME Number of 5-year periods since last cutting 


The first set of screening runs was made on the separate plots (61, 62, 71, 72, and 
Taylor Woods separately) with the purpose of identifying those independent variables most 
highly correlated with the log of basal area growth, and eliminating those independent variables 
that proved to be poor(er) predictors. To do this, the independent variables were classified 
into sets, and then each set was placed in a group (table 24). A set is defined as one or 
more independent variables such that, when chosen by combinatorial rules, all independent 
variables in the set are included (or excluded when not chosen). The only sets containing 
more than one independent variable were those formed with the three independent variables 
LBA. , MBA. , and UBA , , or their squares. 
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Table 24.--Definition of independent variables, groups and sets for the individual plot log 
of basal area growth screening runs 


Variable Group Set : 
Variable 
number number number 

1 1 1 1n(D) 

2 2 1 D 

3 2 Dace 

4 3 pe 

5 4 p2-9 

6 5 D3 

7 6 p39 

8 7 Dt 

9 3 1 In (ARF) 
10 2 In(GRF) 
11 4 1 BAZ 

12 2 BA? 

13 3 BA 

14 4 CCF2 

15 5 CCF 3 

16 6 CCF* 

U7 7 T2 

18 8 T3 
19 9 Tt 
20 5 1 in(PCT) 
“al 2 PCT, 
22 3 PCr 
23 4 1n(DCP) 
24 5 DCP 
25 6 LBA, 
26 MBA, 
27 UBA, 
28 7 (LBA,) 2 
29 (MBA, ) @ 
30 (UBA,) * 
31 8 LBA> 
32 MBA> 
33 UBA> 
34 9 (LBA>) 2 
35 (MBA,) 2 
36 (UBA>) * 
37 10 LBA3 
38 MBA 3 
39 UBA 3 
40 Ha (LBA 3) 2 
41 (MBA 3) 2 
42 (UBA 3) * 
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A group contains one or more sets. The five groups in these screening runs were established 
on the following criteria: group 1 contained the single variable, In(D); group 2 consisted of 
the seven powers on D; group 3 was the two rainfall variables; group 4 was the total stand 
density variables; and group 5 was the position in diameter class sets. Site index was not 
included in these runs because the difference between the site index values on each half plot 
was probably too small to prove significant in modeling. Program REX was then set up to 
tabulate the RMSQR values of all regression equations formed by choosing, at most, one set 
from each group. 


From these screening runs, 12 blackjack and yellow pine equations were selected and the 
regression coefficients determined. These coefficients were next examined for reasonableness 
of behavior and for consistency of performance between plots. In addition, the model with the 
lowest RMSQR was picked for each plot and the residuals of these models were checked for 
normality. To do this, skewness and kurtosis statistics were computed using formulas found in 
Kendall and Stuart (1977). From this examination, I concluded that the residuals are not 
normally distributed, that they are highly skewed, and are leptokurtic (Kendall and Stuart 
1977). One consequence of the lack of normality is that all of the usual significance tests 
that could be made on regression results are not applicable. 


Another problem that made testing difficult or impossible was the lack of data range and 
overlap between the different stand conditions. The site index of the even-aged plots was 
larger than the uneven-aged plots; the stand structures were obviously different; and the even- 
aged plots covered only a narrow range of diameter classes. For the uneven-aged stands, the 
managed and virgin plots differed in stand structure due to cutting. If legitimate tests were 
available, a test indicating significant differences in the models of each stand condition 
might be due to the data being truly incompatible, or merely that the the data are complementary. 
Because of these problems, I concluded that testing to see if individual data sets could be 
pooled was not practical. 


The main question was whether the even-aged and uneven-aged data sets could be combined 
for blackjack pine. To answer this without testing, it was decided to develop two blackjack 
pine equations: one with the even-aged data, and one without. Both of these could then be 
evaluated as to their predictive capability. 


A second set of screening runs was then made with the plots combined as such: for both 
blackjack pine and yellow pine, the virgin, uneven-aged plots (61 and 62) were combined; and 
the cut, uneven-aged plots were combined (71 and 72). In addition, another set of blackjack 
pine data was created by combining the Taylor Woods even-aged data with the cut, uneven-aged 
data. From the findings of the first screening, the independent variables T, PCT, and DCP 
(and their transformations) were eliminated. Site index (S) was added as another group. 


The results of this second set of screening runs were examined, and 11. blackjack pine and 
4 yellow pine equations were selected and the regression coefficients determined. As before, 
these coefficients were also checked for reasonableness of behavior and for consistency of 
performance between data sets. For blackjack pine, those independent variables that both 
behaved reasonably and minimized RMSQR were 1n(D), D, In(S), 1In(GRF), LBA>, MBA>, and UBA>. 
For yellow pine, the independent variables included 1n(D), D, i1n(S), in(GRF), (LBA>)2, (MBA>)2, 
and (UBA,)*. The independent variable of 1n(S) was inconsistent for blackjack pine and 
unreasonable for yellow pine. It was included, however, because of the desire to develop a 
simulator applicable for the range of site indices found in the habitat type. The problems 
with site index persisted throughout the analysis and finally necessitated special handling, 
which will be discussed later. 


The differences in the size of coefficients between the uncut and cut data sets were then 
modeled as a function of time-since-last-cutting (TIME), expressed as the number of 5-year 
intervals since the last cutting. The uncut virgin data sets were given a value of 60 (or 300 
years). I hypothesized that the transition from ''cut'' regression coefficients to "uncut" 
regression coefficients would be along a smooth path, and therefore three sigmoidal curve 
forms were developed using MATCHACURVE techniques (Jensen and Homeyer 1970). All three curves 
are zero when TIME is zero and one when TIME equals 60. The differences in the three curves 
are their forms (fig. 9), and they were chosen more or less subjectively. Because the TIME 
data fell at both extremes, I felt that a more objective approach was not possible. The three 
equations are: 


Ay = -0.244178 + 1.244178*EXP(-(1.176471 - 0.019607*TIME) 2) 

A, = -0.095336 + 1.095336*EXP(-(1.25 - 0.020833*TIME) *) 

Az = -0.000203 + 1.000203*EXP(-(1.428571 - 0.023809*TIME) ©) 
where 


TIME = time-since-last-cutting expressed in number of 5-year growth periods. 
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Time Since Cutting in Five Year Growth Periods 


Figure 9.--Time-since-last-cutting sigmoidal curves. 
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Using these time-since-last-cutting variables, a third set of screening runs was made to 
combine all data sets for each "vigor" class. Again, blackjack pine had two data sets; one 
with Taylor Woods data and one without. Table 25 lists the independent variables and how they 
were grouped. The runs were set up so that the basic independent variables from the previous 
step were forced into all equations. The screening then picked, at most, one independent 
variable from each of the remaining groups. These groups were formed by multiplying the three 
sigmoidal curves (A,, Az, A3) by all the independent variables except 1n(S) and iIn(GRF). 

For those independent variables in which a time-since-last-cutting variable is also picked, the 
effect is to provide for a change to the regression coefficient of that independent variable 
as time-since-last-cutting changes. 


Table 25.--Definition of independent variables, groups and sets for combined plot, log of basal 
area growth screening runs 


Variable Group Set Vomaael 
number number number 
1 il 1n(D) 
2 D 
3 LBAy or (LBAp)? 
4 MBA» or (MBA,)2 
5 UBA» or (UBA) 2 
6 In(S) 
7 1n(GRF) 
8 1 1 A, 1n(D) 
9 2 Ay1n(D) 
10 3 A31n(D) 
11 2 1 AD 
12 2 AoD 
13 3 A3D 
14 3 1 A,LBAy or A,(LBAz) 2 
15 2 AoLBAy or Ay (LBA>) 2 
16 5 A3LBAy or A3(LBAz) 2 
17 4 1 A,MBA> or Aj (MBA) 2 
18 2 A MBAy or A> (MBA)? 
19 3 A3MBAy or A3(MBA,) 2 
20 5 1 A,UBAy or A, (UBA,)* 
21 2 A,UBAy or A> (UBAz) 2 
22 3 A3UBAy or A3(UBA,)* 
23 6 1 Ay 
24 2 Ao 
25 3 A3 


1In those cases where two variables are given, the first is for blackjack pine and second 
is for yellow pine. 


The most promising equations from these runs were selected and the regression coefficients » 
determined. Analysis of these coefficients showed that the coefficients on In(S) were not 
reasonable (that is, they were negative) for yellow pine and for blackjack pine with the 
uneven-aged data set alone (plots 61, 62, 71, and 72), and that the coefficient on in(GRF) was 
not reasonable for blackjack pine with the uneven- and even-aged data sets combined. These 
runs also indicated that a high degree of multicollinearity existed. 
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Multicollinearity and Ridge Regression 


Multicollinearity results when one independent variable is highly correlated with another, 
or with a linear combination of other independent variables. While the model is still an 
unbiased estimator, the effect of multicollinearity is to produce highly imprecise estimates 
of each regression coefficient (Kmenta 1971). This high degree of impreciseness is due to the 
multicollinearity between independent variables causing the correlation matrix to approach 
singularity (Farrer and Glauber 1967). One possible effect is that some of the regression 
coefficients may be of the wrong sign from what was expected (Hoerl and Kennard 1970a, 1970b). 


One method used for reducing multicollinearity effects is ridge regression. This sacrifices 
unbiasness to obtain estimates that, when compared to their unbiased least squares counter- 
parts, are often interpretable and have a smaller mean square error. In terms of solving for 
the standardized regression coefficients,® the method consists of adding a small constant 
value, K, to the diagonal elements of the correlation matrix and then solving in the usual 
manner for the regression coefficients. When K is zero, the ordinary least squares regression 
estimates result. While K can be any positive value, it usually lies between zero and one. 

Also, the larger the value of K, the larger the bias. Details of ridge regression can be 
found in Hoerl and Kennard (1970a, 1970b), Brown and Beattie (1975), Marquardt and Snee (1975), 
or Hocking (1976). 


Numerous methods have been proposed for determining which value of K is "best"; Hocking 
(1976) summarizes these methods. One technique common to many of these methods is the devel- 
opment of a ridge trace, which was first proposed by Hoerl and Kennard (1970a, 1970b). The 
ridge trace is usually a plot of the resulting standardized regression coefficients over a 
range of their respective K values, usually from zero to one. Using the ridge trace, Hoerl 
and Kennard (1970a) suggested four items to consider when deciding upon a value of K: 


1. At a certain value of K, the system will stabilize and have the general characteristics 
of an orthogonal system. 


2. Coefficients will not have unreasonable absolute values with respect to the factors 
for which they represent rates of changes. 


3. Coefficients with apparently incorrect signs at K = 0 will have changed to have the 
proper sign. 


4. The residual sum of squares will not have been inflated to an unreasonable value. It 
will not be large relative to the minimum residual sum of squares or large relative to what 
would be a reasonable variance for the process generating the data. 


Because of the high degree of multicollinearity in the data and of the unreasonable signs 
on some of the independent variables, I decided to try ridge regression to determine if it 
could improve the estimates of the regression coefficients. I wrote a program to compute 
regular and standardized regression coefficients for K values ranging from zero to 0.95 in 
intervals of 0.05. These values could then be plotted to produce the ridge trace. 


6Standarized regression in coefficients, b.*, and normal regression coefficients, b., are 
related in the following manner: 
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The results for yellow pine indicated that most of the independent variables were stable. 
Two exceptions are the independent variables D and In(D). Their behavior was common to all 
data sets and to all models, and was expected because of the high correlation between the two. 
Disappointingly, ridge regression did not change the sign on In(S). For blackjack pine fitted 
to the uneven-aged data alone, the results were similar. The sign of In(S) did not change, 
and the four independent variables involving D all showed instability. This last was the 
most notable feature of the ridge trace for blackjack pine fitted to both the even- and uneven- 
aged data sets. 


Because of the desire to have a model that was applicable over the range of site indices 
found in the habitat type, the results of the ridge regression runs were unacceptable. I 
decided, therefore, to force a reasonable value of In(S) on the models, which was done by 
fitting a new dependent variable of In(BAG/S). This has the same effect as forcing the 
independent In(S) onto the log of basal area growth model with a regression coefficient of 
one. The choice of one as the regression coefficient was made because, fixing other factors, 
growth increased in direct proportion to an increase in site index (W. H. Meyer 1938; Larson 
and Minor 1968; Ek 1974; and Myers and others 1976). 


Given this new dependent variable, a final set of screening runs was made utilizing the 
same strategy as that used for the third set. Promising models were chosen from these runs 
and additional runs were made to determine their regression coefficients. Ridge regression 
was then used to examine the stability of the various independent variables. 


As with the ridge trace for the log of basal area growth, the independent variables 
showing the greatest amount of instability are those involving D. For both blackjack pine 
equations, a moderate amount of instability was also exhibited by those variables involving 
the same basic variables (for example A3*UBA, and UBA)). 


The instability between D and In(D) was expected because of their high correlation, but 
this was acceptable because the two independent variables are necessary to provide the desired 
"neaking'' effect over diameter class size. For those independent variables involving time- 
since-last-cutting, however, I decided to eliminate those that were highly unstable and 
therefore did not greatly improve the least squares regression equations. This resulted in no 
independent variables being removed for yellow pine, but it did result in three independent 
variables being removed from both blackjack pine equations. The ridge traces for the finished 
equations are found in figures 10, 11, and 12. 


For two reasons, the independent variable In(GRF) was also eliminated. First, the fact 
that the variable was negative for the blackjack pine equation with the even-aged data was 
cause for concern. The even-aged data were collected 40 years later than the uneven-aged 
data. The sign change caused by the incorporation of the even-aged data might therefore 
signal that the influence of rainfall was not uniform over time. Also, the data needed to 
compute In(GRF) are not readily available. The data source used in this study came from the 
Fort Valley Experimental Forest, which has a weather station. 
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Figure 10.--Final ridge trace (with growth period rainfall) for yellow pine using the 
dependent variable of In(basal area growth/site index). 
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Figure 11.--Final ridge trace (with growth period rainfall) for blackjack pine fitted 
to the uneven-aged data alone with the dependent variable of In(basal area growth/ 
site index). 
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Figure 12.--Final ridge trace for blackjack pine fitted to the even- and uneven-aged 
data with the dependent variable of In(basal area growth/site index). 
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APPENDIX D 
CORRECTION FOR LOG BIAS 


If the assumption of normality of residuals cannot be met, then correction of a log model 
for bias can be a problem for which no easy solution is apparent. The following, however, is 
a proposed correction factor based more on intuition than mathematical rigor: 


If S is such that 1nS\N(u,o2) then it can be shown that: 
E[S] = EXP[u + 40%] (Mood and others 1974) 
Therefore, an estimator of S is: 
S = EXP[in(S) + MSE] 
This correction for lognormal bias, 4:*MSE, is the same as that proposed by Oldham (1965) and 
Baskerville (1972) as an approximate correction to the log regression model. Now, if S is not 
lognormally distributed, it would still seem reasonable to assume that a constant, K, still 


exists such that: 


S = EXP[in(S) + K] 


Rearranging, 
S = EXP[in(S)] EXP[K] 
EXP[K] = : 
EXP[1n(S) ] 


i= ia) = = ale 
The regression equation: 

DY) = <bgete bye) + 2c + bx, 
can also be expressed as: 

In(¥) = In(Y = by (x, = x1) + Bo(xo = xo) + 22: + De Gae <) 
or, 


In(Y) = iIn(Y) + by (x = x) cf bo (xo = xo) ap OOO ar b (x, x) 
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If K = I1n(Y) - In(Y) were added to this model as an approximate correction for log-bias in 
the nonnormal case, then the effect would be: 


In(¥) = K + In(¥Y) + by (x, - x) + bo(x2 = x2) fa rece b (x, if x,) 


& In(Y) = In(Y) + ini(y) + by (x, - x1) EE heey a bx, - =) 
pic 1n(Y) iff by, (x} is x1) + bo (x9 - x») + see + boy, = x) 


In other words, the correction would force the model through the log of mean Y instead of the 


mean of log Y. Therefore, Y will pass through Y when x; = x1, x2 = x2, °**, 5 which 


appears to be the adjustment desired. Whether this is a true 'mean-unbiased" estimator is 
impossible to prove without first knowing the true distribution of In(Y). 


The correction for log bias is applied by adding it to the intercept term of the log of 
basal area growth model. The log bias correction factors for the three equations developed 
in this study are found in table 11. For comparative purposes, the appropriate correction 
for normally distributed residuals (that is, %-MSE) is also found in table 26. Notice that 
if normality had been assumed, the correction would have been considerably larger than the one 
not assuming normality. 


Table 26.--Proposed correction factor and log-normal correction factor for log bias of the 
three log of basal area growth equations 


Equation Proposed factor Lognormal factor 
type [ln(mean Y)] - [mean In(Y)] 15>MSE 
Blackjack pine with 
even-aged data 0.300534069 0.48453144 
Blackjack pine without 
even-aged data - 324720185 -56464436 
Yellow pine - 873936593 2.76746302 
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APPENDIX E 


DEVELOPMENT OF MORTALITY MODELS 
Prior Findings 


For uneven-aged stands, noncatastrophic mortality has been found to be correlated to 
diameter size, severity of cutting, and whether the trees were blackjack or yellow pines. 
Krauch (1926) reported that percent mortality increased as diameter increased. When percent 
mortality is plotted over diameter, the curve is U-shaped with the minimum near 20 inches 
(Pearson 1939, 1950; Pearson and Wads:orth 1941; Wadsworth and Pearson 1943; and Myers and 
Martin 1963). 


For severity of cutting, Krauch (1926) discovered that "the percentage of mortality 
increases in inverse ratio to the degree of cutting.'' This discovery was supported by Lexen 
(1935). Wadsworth and Pearson (1943) found that percent mortality of blackjack or yellow pine 
alone was higher in a cutover stand than in a virgin stand, but that this situation reversed 
itself when the two were combined. They attributed this reversal to peculiarities of percent- 
ages. Pearson (1950) also found that a virgin stand had a higher rate of mortality than a 
moderately cut stand, but that severe cutting increased mortality because of windthrow and 
lightning losses in the few large trees remaining. Higher mortality rates in yellow pine than 
in blackjack pine have been reported by Pearson and Wadsworth (1941) and by Wadsworth and 
Pearson (1943). 


Mortality should also be related to productivity, total stand density, and time-since-last- 
cutting. For even-aged stands, Myers and others (1976) found that site and total stand basal 
area were significant in predicting total stand mortality, and there is no reason to believe 
this effect is not also significant in uneven-aged stands. 


While it seems intuitive that the longer the time-since-last-cutting, the higher the 
mortality rate, both Lexen (1935) and Pearson (1939, 1950) have reported that they could find 
no consistent relationship between mortality and time-since-last-cutting. Their findings 
could have been clouded by other effects or by the fact that the length of time-since-last- 
cutting they examined (20 years or less) was not long enough for the effect to manifest itself. 


Definition of Variables and Transformations 


A mortality rate can be expressed in several ways. Ek (1974) and Adams and Ek (1974) 
expressed it directly as the number of trees dying in a diameter class during the growth 
period. Lee (1971) expressed mortality as a percentage of the total stand dying per year. 
Finally, Hamilton (1974), Hamilton and Edwards (1976), and Monserud (1976) all expressed 
mortality as the probability (or proportion) of a tree dying in a growth period. 


The form of the equations has also differed. Ek (1974) used nonlinear least squares 
techniques, while Lee (1971) and Adams and Ek (1974) used ordinary, linear least squares 
techniques. Hamilton (1974), Hamilton and Edwards (1976), and Monserud (1976) all used the 
logistic function fitted using weighted, nonlinear least squares techniques. 


Expressing mortality as a percentage or proportion has an advantage over expressing it as 
the number of trees dying. If properly modeled, a proportion (or percentage) is bounded by 0 
and 1 (or 0 and 100). Multiplying this rate by the number of trees in the class will always 
give a value less than or equal to the number in the class. As a result, predicted number of 
survivor trees is never negative; a result possible under the other approach. 
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A problem can arise with using ordinary, linear least squares regression for predicting 
the percent or proportion of trees dying. Linear models can cause the predicted values to 
exceed zero or one, especially if the prediction falls outside the range of values in the 
developmental data set (Hamilton 1974; Hamilton and Edwards 1976). However, the nonlinear 
logistic function, 


y = (120° paleo RPA asa ao aire 2 ay 


limits the value of y to between zero and one. Also, when the dependent variable is dichotomous 
(which mortality is), the logistic function appears to have improved statistical properties 

over ordinary, linear least squares regression models (Hamilton 1974; Hamilton and Edwards 
L976) 


Hamilton (1974) presents a computer program, RISK, for determining the parameters of the 
logistic function using weighted, nonlinear least squares regression techniques. In using 
RISK, the dependent variable is assigned a value of zero if the tree survived to the end of 
the growth period, and a value of one if it did not. By using only diameter class or stand 
attributes as independent variables, the result is the probability of a tree dying under the 
specified diameter class and stand conditions. 


The reported U-shaped curve of percent mortality over diameter was for all ponderosa 
pine. By separating ponderosa pine into the two vigor classes of blackjack pine and yellow 
pine, the trends will differ from the U-shaped curve. For blackjack pine, it is possible that 
percent mortality will decline as diameter increases; while for yellow pine, the opposite is 
expected. Because blackjack pine usually represents the small, vigorous trees and yellow pine 
the large, overmature trees, it is easy to see why the resulting combined graph could be 
U-shaped. 


Any transform on diameter class size (D) that causes the transformed value to increase as 
D increases would allow the desired effect in the logistic function. The two independent 
variables chosen, D and D*, both meet this criterion and are also simple and interpretable 
(D is straightforward, and Dowes analogous to’ diameter class basal area). As with the basal 
area growth equations, vigor was incorporated into the equations by developing separate 
mortality equations for blackjack pine and for yellow pine. 


Two measures of productivity were tried: site index for the one-half plot (S), and average 
5-year growth period rainfall (GRF). Mortality will decline as productivity increases; 
therefore, the independent variables chosen were the untransformed values S and GRF. A 
positive sign on the regression coefficients of these independent variables would provide the 
desired om j,expected effect; 


The choice of the independent variables to represent stand density was based, in part, on 
the results of the growth analysis. The sets of LBAz, MBAz, UBAz, and of LBA»*, MBA», UBAd?, 
proved best in the growth analysis at representing both the effect of stand density and of 
diameter class position. Also chosen for analysis in this phase were the independent variables 
of total stand basal area (BA) and percentile in the basal area distribution (PCT). 


For time-since-last-cutting, two of the three transforms used in the growth analysis (A, 
and A3) were also tried in this analysis. Expressing severity of cutting was a problem. On 
the plots with cutting, the recording of data.did mot start until) atter jtheycuc. whnuse 
information concerning the amount of trees removed was not available as an index of severity 
of cutting. While the amount of the residual stand might provide some insight into the severity 
of cutting, it is not the best measure. Because stand density reflects the amount of the 
residual stand, it was felt that additional measures were neither justified nor necessary. 
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With the inherent difficulty of predicting mortality, I decided to include an independent 
variable that incorporated all of the above independent variables into one. The obvious 
choice for independent variable is predicted basal area growth. The probability of mortality 
would decrease as basal area growth increased. Predicted basal area growth was used instead 
of actual basal area growth because all independent variables must be observable when the 
model is actually applied in the stand simulator (Monserud 1976). The use of actual basal 
area growth in equation development and then predicted basal area growth in application of 
the equation would introduce a random variable as an independent variable, which violates one 
of the assumptions of least squares regression. To clarify this last point, consider the 
following example: 


Suppose you wanted to obtain an estimator of Y as a function of X (the predicted value 
of X) because once the estimator of Y is determined, the only value available to activate the 


estimator is X. However, the observed values (Y, X) do exist prior to the development of the 


estimator. Now if Y and X are related by 


Yo= Bo. Bik 40<) (1) 
and if 

ou Eo Or (2) 

ee (3) 


then, if X (the observed value) is used instead of X (the predicted value) to develop equation 
foipetheceLtect ss to model 


x¢ Bo + B, (X-e€2) ur eS] (3) into: (@) 


Bo * ByX-Byeo + ©] 


and Kmenta (1971) has shown that the effect of this is to make the estimators of By and 8) 
inconsistent. 


The three predicted basal area growths used were: blackjack pine with the even-aged data 
(BJBAG1), blackjack pine without the even-aged data (BJBAG2), and yellow pine (YPBAG). 


“Screening” Independent Variables 


Program RISK helped determine which of the independent variables best predicted mortality. 
To do this, a run was created to fit 14 equations for the three blackjack pine data sets and 
11 equations for the two yellow pine data sets. The three blackjack pine data sets were (1) 
virgin uneven-aged data, (2) managed uneven-aged data, and (3) managed even-aged data; and the 
two yellow pine data sets were (1) virgin uneven-aged data, and (2) managed uneven-aged data. 
These preliminary data sets were chosen to both reduce the expense of running RISK on large 
data sets, and to allow an examination and subsequent modeling of the trends in model param- 
eters across the various data sets as functions of stand structure and time-since-last-cutting, 
if necessary. 


For both blackjack pine and yellow pine, the results of this "screening" run indicated 
that the only independent variables with predictive strength were D, D?, and the basal area 
growth equations. Because the variable D2 proved to be better than D, it was used in 
the remainder of the analysis. 
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The next step was modeling the change in model parameters across time-since-last-cutting 
and stand structure. To model the effect of stand structure, two new variables were intro- 
duced: RMBA, and RUBA,, which were computed by dividing MBA, (or UBA.) by total stand basal 
area. Thesé variableS are based on the hypothesis that an éven-aged stand would exhibit more 
basal area in the middle basal area class (or less of its total basal area in the upper basal 
area class) than would an uneven-aged stand. 


Utilizing the time-since-last-cutting and two of the "stand structure" variables (RUBA3 
and RMBA>) as independent variables themselves or as multipliers with the other variables, a 
second set of ''screening'' runs was made on three data sets: blackjack pine with the even-aged 
data, blackjack pine without the even-aged data, and yellow pine. From these runs, the two 
most promising models from each data set were selected for more detailed analysis. The 
independent variables for these models are found in table 27. Both blackjack pine growth 
equations were used in each of the two blackjack pine data sets because no choice has been 
made as to which blackjack pine growth equations will be final. 


Table 27.--Independent variables in mortality models selected from each data set for final 


analysis 
Data set Model number Independent variables 

Blackjack pine with 1 BJBAG1, RMBA> 
even-aged data 

2 BJBAGI1, RUBA3 

3 BJBAG2, RMBA> 

4 BJBAG2, RUBA3 
Blackjack pine without 1 BJBAG], A»*BJBAG1 
even-aged data 

2 A,*BJBAG1, D2 

3 BJBAG2, A, *BJBAG2 

4 A, *BJBAG2, D? 
Yellow pine 1 YPBAG, A, *YPBAG 

2 YPBAG, A,*YPBAG, D2 


Up to this point, the selection of the ''best'' model has been made by jointly considering 
four factors: (1) reasonableness of signs on the independent variables; (2) overall signifi- 
cance of the model as indicated by an "'F" statistic in RISK; (3) maximization of "'t" values on 
each parameter, as computed by RISK; and (4) minimization of a chi-square ''goodness-of-fit" 
value also computed by RISK. Of the three measures of fit computed in RISK (R, t, and chi- 
square), Hamilton (1974) preferred the chi-square. Unfortunately, the chi-square tables in 
RISK used predicted event classes too wide for the mortality data used in this study. Asa 
result, the observations fell in very few of the classes. To correct for this, a program was 
written to compute the chi-square statistic over a reasonable range of predicted mortality 
classes and also over diameter classes. 


Using this program, chi-square values were determined for the models selected for final 
analysis. Analysis of the results indicated that both yellow pine equations appear to perform 
well. Because model 2 (the model with YPBAG, A,*YPBAG and D2) has slightly better chi-square 
values, it was chosen as the best yellow pine mortality model. The goodness-of-fit test 
across diameter classes is found in table 10. 
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The results for the blackjack pine equations developed using just the uneven-aged data 
were also satisfactory. The chi-square values for all the tests were not significant at the 
99 percent level (the testing level). The signs on the coefficients of models 1 and 3 were, 
however, of opposite values. This unexpected result caused the following: for stands recently 
cut, the higher the basal area growth rate, the higher the mortality rate; while for virgin 
stands, the opposite was true. The coefficients on models 2 and 4 were reasonable, and so 
they were picked as the best models for uneven-aged blackjack pine. Tables 11 and 12 provide 
the goodness-of-fit tests across diameter classes for the final models. 


The chi-square values indicated poor fits for the equations developed using both the 
even- and uneven-aged data, and this was disappointing. So a new group of runs was tried. 
The independent variables used in these runs are found in table 28. Of these, the two best 
sets of equations were models 19, 20, and models 21, 22. While the models are an improvement, 
they are still significantly different at the 99 percent testing level used. The goodness-of- 
fit tests across diameter classes for the best models (models 21 and 22) are found in tables 
13 and 14. 


Table 28.--Independent variables used in final set of runs for even- and uneven-aged blackjack 
pine mortality 


ae Independent variables 
5 BJBAG1, D*, A,, LBA>, MBA>, UBA> 
6 BJBAG2, D2, A,, LBAz, MBA>, UBA> 
7 BJBAG1, D2, A,, MBA» 
8 BJBAG2, D2, A,, MBA» 
9 BJBAG1, A,, MBA» 
10 BJBAG2, A,, MBA», 
at BJBAG1, D2, A,, A,2, MBA, 
12 BJBAG2, D2, A,, A,2, MBA» 
13 BJBAG1, D*, TIME, TIME*, MBA, 
14 BJBAG2, D*, TIME, TIME*, MBA, 
15 BJBAG1, A,, A;*, MBA» 
16 BJBAG2, A,, A,*, MBA» 
17 BJBAG1, TIME, TIME*, MBA, 
18 BJBAG2, TIME, TIME2, MBA 
19 BJBAG1, D2, TIME, TIME, MBA> 
20 BJBAG2, D2, TIME, TIME’, MBA> 
21 BJBAG1, D2, TIME, TIME -, MBA> 
22 BJBAG2, D2, TIME, TIME~~, MBA> 
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Examining causes for mortality provides insight into the problem of modeling the combined 
uneven- and even-aged mortality data. For the uneven-aged data, the principal causes of 
mortality were lightning, wind, bark beetles, and mistletoe (Pearson 1939; Wadsworth and 
Pearson 1943; Myers and Martin 1963). The amount would indicate that these are noncatastrophic 
forms of mortality. The main causes for mortality in the even-aged plots were snowbreak and 
rust (Schubert 1971, 1974). Of the 160 trees that died on those subplots used in equation de- 
velopment, 86 percent died because of snowbreak in the first growth period after cutting, and 
12 percent died of rust in the second period. 


A problem exists as to whether snowbreak should be considered noncatastrophic or catastrophic 
mortality. Accelerated mortality after cutting is not unique to the Taylor Woods study. Also, 
windthrow after cutting has been widely reported (Alexander 1973 and 1975; Reukema 1970; 

Reukema and Pienaar 1973; and Williamson and Price 1971) while both Boldt (1970) and Williamson 
(1976) reported serious problems with snowbreak. Barrett (1965) also noted that "extensive 
precautions were taken to protect the stand from damage by porcupines and snowbend, the two 
greatest threats to maintaining the designated tree densities.'' The problem of snowbreak and 
bending has been described as being common to the Southwest, particularly among young dense 
stands such as the Taylor Woods subplots (Pearson 1950 and Schubert 1974). 


Should, therefore, a mortality agent described as common be considered as a natural 
result of thinning? The answer is subjective and cannot be presented satisfactorily in this 
study. To examine how the inclusion of the Taylor Woods mortality data affects the prediction 
of uneven-aged mortality, the chi-square goodness-of-fit tests were made to the even- and 
uneven-aged mortality equations using just the uneven-aged data. The resulting chi-square 
statistics indicated that the inclusion of the even-aged data has deteriorated the predictive 
capability for the uneven-aged plots. Whether this result indicates that snowbreak is non- 
normal (that is, catastrophic) mortality is not clear. A comparison of the goodness-of-fit 
tests across diameter class (tables 29 and 30) to those for the uneven-aged equations (tables 
11 and 12) indicates that the greatest differences occur in the 4- through 6-inch diameter 
classes, and these are the classes where most of the Taylor Woods data occurred. Despite 
these problems, I decided to also test these two equations in the validation phase along with 
the uneven-aged blackjack and yellow pine equations. Perhaps the equations developed with the 
even-aged data will prove to be better predictors over the long run. 
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Table 29.--Chi-square test across diameter classes for best even- and uneven-aged blackjack 
pine mortality equation (using basal area growth equation for even- and uneven- 
aged blackjack pine) compared to uneven-aged data only 


Diameter Number of trees Actual Predicted Chi-square 


class in class mortality mortality value 
4 4,416 11 45.57 26.2 
5 4,066 18 35.90 8.9 
6 Seoze 18 24.72 138 
7 2,629 iy 17-62 0 
8 2,021 8 11.44 TO 
9 15551 7 7.84 Cer 
10 1,203 5 5ED8 yf 
11 1,004 6 4.22 ah 
12 893 3 3.54 api 
13 737 2 IDS TIS «2 
14 686 1 22535 .9 
15 666 6 2.61 4.4 
16 685 4 2265 oa, 
U7, 649 2 2252 cal 
18 610 1 2.38 .8 
19 473 0 1.96 2.0 
20 - 21 908 6 4.07 39 
22 =925 622 9 3.06 LSS 
24 - 25 553 1 2213 .6 
26+ 261 0 2250 25 


wn 


Chi-square = 6 


Ws 


Table 30.--Chi-square test across diameter classes for best even- and uneven-aged blackjack 
pine mortality equation (using basal area growth equation for uneven-aged black- 


Diameter 
class 


20) = 2:1 


22 = 23 


24 = 25 


26+ 


4 


5 


6 


jack pine only) compared to uneven-aged data only 


Number of trees 
in class 


4,416 
4,066 
3,322 
2.629 
2021 
1,551 
1,203 
1,004 


893 


$55 


261 


Actual 
mortality 


Wi 
18 
18 


7, 


74 


Predicted 
mortality 


44, 


Sys 


24 


17. 


It 


93 


80 


Kisii 


74 


55 


Bo 


.65 


5249) 


.62 


-82 


.61 


aA 


is 


5 (015) 


-49 


.05 


215 


. 03 


. 04 


. 34 


Chi-square 
value 


ADSEO 


8.9 


1Vs8 


Chi-square = 62.8 


APPENDIX F 
DEVELOPMENT OF A CONVERSION MODEL 


I hypothesized that the following independent variables might prove useful for predicting 
the conversion rate: diameter class size (D), site index (S), total number of trees 7.6 inches 
and larger (TT8), total basal area in trees 7.6 inches and larger (BA8), total number of trees 
in the diameter class (TTDC), and total basal area in the diameter class (BADC). Except for 
site index, I expected that as each of these independent variables increased, the conversion 
rate would increase. The opposite would hold true for site index. 


The expected relationships between the proposed independent variables and the conversion 
rate were based upon how changes in these independent variables would affect growth. Growth 
and vigor are intimately related, and it was expected that a decrease in growth rate would 
cause an increase in the conversion rate. Therefore, any of the proposed independent variables 
that would reduce the growth rate as they increased in value would also be expected to increase 
the conversion rate. 


Based on the same reasoning as discussed in the mortality rate section, I decided to 
express the conversion rate as a proportion. To model this, program RISK was again used. The 
dichotomous dependent variable was defined as 1.0 if the tree converted to a yellow pine, and 
a value of 0.0 if it did not. 


Twenty-two regression runs were then made using the data reserved for modeling. These 
runs were selected to cover most, if not all, reasonable combinations of the independent 
variables. The runs allowed decisionmaking at various points during the analysis on what 
combinations of independent variables to examine next. Chi-square values across diameter 
classes and across predicted conversion classes for the three best models were then examined. 
In all cases, the predicted values were significantly different from the actual values at the 
99 percent testing level. An analysis of the chi-square values for the best model indicated 
the primary area of misfit was in the 12-inch diameter class. A review of the basic data also 
revealed that all of the conversions in the 11-, 12-, and 13-inch diameter classes and two of 
the three conversions in the 14-inch diameter class occurred on subplot 13 of plot 61. 


I then decided to circumvent the problem of lack of fit by strengthening the data base. 
This was done by adding to the conversion data base those subplots originally eliminated 
because of a highway running through them. Using this expanded data set, the 25 new runs were 
made using the same procedures as described for the first set of runs. Again, the chi-square 
values for the three best models indicated a significant difference between predicted and 
actual values. As before, the lack of fit was due to the high number of conversions on subplot 
13 in the 11-, 12-, 13-, and 14-inch diameter classes. 


I eliminated subplot 13 from the data base as being atypical of the conversion process, 
then made the 10 runs listed in table 31. This time, the chi-square values of the three best 
models (5, 9, and 10) indicated no significant difference at the 99 percent testing level 
between predicted and actual conversion values across predicted conversion classes. Across 
diameter classes, only model 9 did not differ significantly from the actual data (table 16) 
and was therefore chosen as the final model fitted to this data set. 


To test model 9 further, the model was checked against validation data. Across predicted 
conversion classes, the chi-square value for the validation data was 140.9 with 6 d.f., while 
across diameter classes it was 178.8 with 7 d.f. In both cases, the difference between predicted 
and actual conversion was significant at the 99 percent testing level. An examination of the 
validation data showed that, like subplot 13 of the original data, subplot 16 exhibited a 
large number of conversions in the 12-inch diameter class. Why this peculiarity exists in 
this diameter class on these two subplots is unknown, but I eliminated subplot 16 from the 
validation data also. The resulting chi-square fit for this modified validation data set was 
3.4 with 7 d.f. across diameter classes (table 32) and 3.4 with 6 d.f. across predicted con- 
version classes. Both of these chi-square statistics were insignificant, indicating a good 
fit to the modified validation data set. 
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Table 31.--Listing of conversion models tested using expanded data set minus subplot 13 


ae Independent variables 
1 DESH SHE BAS 
2 Det Ses BADE 
3 Dy, Say LEDC. BAS 
4 D, S, TTDC, BADC 
5 De G08 
6 DA uEDG 
7 D2, TTS 
8 D2, TTDC 
9 D, D2, S 
10 D, S 


Table 32.--Validation chi-square test across diameter classes for final conversion model 


Diameter Number of trees Actual Predicted 

class in class conversion conversion 
4-15 747 2.00 eZ 
16-19 84 1.00 Ibs oy! 
20-21 42 1.00 ISX) 
22 24 1.00 OS 
23 26 1.00 19358 
24-25 23 1.00 Weed 
26-28 13 3.00 28) 
29+ 11 2.00 L365 


Chi-square statistic 
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APPENDIX G 


DEVELOPMENT OF RECRUITMENT MODELS 
Prior Findings 


Natural regeneration is the process in which an adequate supply of viable seed is produced, 
dispersed over the area, germinated, and then survives to produce a full stocking of established 
seedlings. The amount of seed a tree produces depends upon tree size, vigor, level of competi- 
tive stress, and presence of cone-damaging agents (Larson and Schubert 1970; Schubert 1974) as 
well as environmental factors. Schubert (1974) concluded that good cone crop years occur at 
intervals of 3 to 4 years. The amount of viable seed is positively correlated with size of 
cone crop (Larson and Schubert 1970; Schubert 1974). 


Seed is disseminated in the fall and generation occurs the following summer (Schubert 
1974). Germination follows adequate rainfall and, if this takes place too late in the season 
or fails to take place at all, mortality rates can be extremely high (Pearson 1950; Schubert 
1974). Second-year seedlings are also highly susceptible to drought. Therefore, successful 
regeneration requires adequate rainfall in two consecutive years (Meagher 1950). Seedling 
survival is also heavily influenced by the amounts of competing vegetation (Pearson 1942, 
1950; Schubert 1974). 


The ingrowth rate, therefore, should be influenced by the potential of the stand to 
produce cones, which is positively correlated with tree size (Larson and Schubert 1970), by 
the level of competition, and by the time since the last good seedling year. Another factor 
that might influence the ingrowth rate is the structure of the competition. For example, an 
even-aged stand and an uneven-aged stand could have the same potential for producing cones and 
the same overall level of competition, but because the diameter structures of the two stands 
are different, the ingrowth rates might differ. (The same total basal area spread over many 
diameter classes in the uneven-aged stand may produce a different level of competition on a 
given diameter class from the same total basal area concentrated in a narrow range of diameter 
classes in an even-aged stand.) 


The only whole-stand ingrowth models found in the literature were the two tried by Moser 
(1972, 1974) and one reported by Ek (1974). Moser's differential equations were of the form: 


=P2QMD 


d(Ingrowth)/dt = b, (1) (Moser 1972) 


and 


bo + ee (2) (Moser 1974) 


d(Ingrowth) /dt 
where 
QMD = quadratic mean stand diameter 


BA total stand basal area. 


Integrating these functions over the time interval tg to (tp + 5) provides a means for con- 
verting these differential forms to equations that represent ingrowth in any 5-year growth 
period. The integrated equations are: 


Th 


Ingrowth = t+*bj,e 


b>QMD =b2QMD 


- tob) 


(3) 
and 


Ingrowth = bot + b,BA 


= 5bo + 5b,BA 
= be + D*BA (4) 


Ek's (1974) model predicts the total number of ingrowth trees directly by the equation: 
by .~-1 
Ingrowth = Bale. re 


(5) 


where 


alt 


I 


total number of trees 


BAY ="total’ sittand ‘basal area. 


Moser used nonlinear regression to fit model (1) and linear regression to fit model (2). 
Ek also used nonlinear regression to fit model (5). Both are assuming, therefore, that the 
error term is additive, but neither provided evidence to support this assumption. 


Total Ingrowth Model Development 


The previously discussed literature review for ponderosa pine in the Southwest indicated 
that ingrowth might be influenced by the potential of the stand to produce cones, the level 
and structure of competition, and the time since the last good seedling year. 


From the data presented by Larson and Schubert (1970), the following equation was devel- 
oped to predict the number of cones a tree in a given diameter class (greater than or equal to 
12 inches) would produce in a year: 

Number of cones per tree per year = 0.94993974(D-12) + 0.61427282(D-12) 2 

Deel 
where 

D = diameter class size. 

Using this equation, the predicted number of cones (C) for the stand was computed and 
used as an independent variable. The number of cones is expected to be positively 
correlated with ingrowth. 

To represent level and structure of competition, six diameter classes were created (4 to 
6 inches, 7 to 9 inches, 10 to 1S’inches, 16 to) 21) inches, 22 to: 27 inches.) 28-4 sinches)), and 


the number of trees and basal area in each was determined. These classes can then be combined 
in numerous ways to produce various measures of the level and structure of competition. 
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Three additional diameter classes were also created by determining the 1l-inch diameter 
class in which the mean stand diameter exists and then adding to it the surrounding one, two, 
or three l-inch diameter classes. The number of trees and the basal area in each of these 
three "mean stand diameter classes" were then determined. 


Because of insufficient data, the time since the last good seedling year was not incor- 
porated into this study. 


Other independent variables thought to be important for predicting ingrowth were site 
index, mean stand diameter, and quadratic mean stand diameter. Moser's (1972) work suggested 
the latter two independent variables. The definitions of all variables used to predict total 
ingrowth are found in table 33. 


In the first phase of the regression analysis, a set of screening runs was made on the 
virgin, uneven-aged data and on the managed, uneven-aged data using program REX. The variables 
were divided into sets and groups as defined in table 34, and the combinations selected were 
made by picking, at most, one set from each group. This resulted in 6,334 regressions being 
examined for each data set. 


From these screening runs, seven of the most promising models were selected, their 
regression coefficients were determined, and then examined for reasonableness of behavior 
between data sets. A preliminary analysis was also made concerning normality of residuals. 
The skewness and kurtosis statistics indicated that the residuals were not badly "nonnormal," 
and therefore "t'' tests could be used to roughly check the significance of the regression 
coefficients. 


Based on the reasonableness and significance criterion, the six original basal area 
eiasses (BC, , BE, BC.) BC,,. BE, BC.) were collapsed, first to the three classes (BC,; BCs 
+ BC3; BC, + BCs + BCg) and ultimately to two classes (BC,; BCy + BC3 + BC, + BCs + BCg). 
The latter two classes were redefined as BACL, and BACL>. 


The next step was to combine the virgin and managed data sets. The two best models from 
the previous set of runs were chosen, and two new screening runs were made. These runs "forced" 
the basic models upon the combined data sets and then allowed adjustment for time-since-last- 
cutting by screening all combinations of the products of the independent variables with the 
three time-since-last-cutting variables (A;, Aj, and A3). (The definitions of groups, sets, 
and variables are found in table 35.) 


An examination of the runs disclosed that the signs of the coefficients on both In(site 
index) and 1n(4-inch diameter class basal area growth) were negative. This was not considered 
reasonable. The desire to have site index as a component of the total ingrowth model dictated 
that a reasonable function of site index be "forced"' upon the model by fitting 1n(Total 
Ingrowth/Site Index) as the dependent variable. This approach was the same as used earlier in 
developing the basal area growth equations. Another screening run was performed using the 
Same strategy as for model 1 of table 35; except the In(site index) independent variable was 
eliminated. From this run, I concluded that the basal area classes could be further collapsed 
to the final two classes BACL, and BACL». A final screening run was then made to provide the 
following log model: 


m2 


Table 33.--Definition of variables used to model ingrowth 
Variable Definition 


BJBAG1-4" Predicted basal area growth of the 4-inch diameter class 
using the even- and uneven-aged blackjack pine equation 


BJBAG2-4" Predicted basal area growth of the 4-inch diameter class 
using the uneven-aged blackjack pine equation 


S Minor's (1964) site index 

C Predicted number of cones per acre per year 

MD Mean stand diameter 

QMD Quadratic mean stand diameter 

TIME Time-since-last-cutting 

ly Total number of trees per acre 

BA Total basal area per acre 

BC, Basal area in the 4-inch through 6-inch diameter classes 

BCo Basal area in the 7-inch through 9-inch diameter classes 

BC3 Basal area in the 10-inch through 15-inch diameter classes 

BCy Basal area in the 16-inch through 21-inch diameter classes 

BCs Basal area in the 22-inch through 27-inch diameter classes 

BCg Basal area in the 28+-inch diameter classes 

Ty, Total number of trees in the 4-inch through 6-inch diameter classes 
TC» Total number of trees in the 7-inch through 9-inch diameter classes 
TC3 Total number of trees in the 10-inch through 15-inch diameter classes 
TCy Total number of trees in the 16-inch through 21-inch diameter classes 
TGs Total number of trees in the 22-inch through 27-inch diameter classes 
T¢ Total number of trees in the 28+-inch diameter classes 

BMD, Basal area in the diameter class of mean stand diameter plus the 


adjacent diameter classes 


BMD > Basal area in the diameter class of mean stand diameter plus the 
adjoining four diameter classes 


BMD 3 Basal area in the diameter class of mean stand diameter plus the 
adjoining six diameter classes 


TMD, Number of trees in the diameter class of mean stand diameter plus 
the adjacent diameter classes 


TMD» Number of trees in the diameter class of mean stand diameter plus 
the adjoining four diameter classes 


TMD 3 Number of trees in the diameter class of mean stand diameter plus the 
adjoining six diameter classes 
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Table 34.--Specification of groups, sets, and variables for first set of total ingrowth 
Screening runs 


Variable Group Bee Variable 
number number number 

1 1 1 In(BJBAG1-4') or 

1n(BJBAG2-4") 

2 2 1 In(C) 

3 3 il 1n(S) 

4 4 l 1n(MD) 

5 2 1n(QMD) 

6 5 1 In(T) 

7 2 1n(BA) 

8 3 In(BC,) 
9 In(BCo) 
10 1n(BC3) 
el 1n(BC,) 
12 1n(BCs) 
15 In(BCg¢) 
14 4 1n(BMD, ) 
15 5 In(BMD >) 
16 6 1n(BMD3) 
17 7 In(TCq) 
18 In(TC>) 
19 In(TC3) 
20 In(TCy) 
21 In(TCs) 
22 1In(TCg) 
23 8 1n(TMD,) 
24 9 1n(TMD>) 
25 10 1n(TMD3) 

26 6 1 1p 

27 2 BA 
28 3 BAl-° 
29 4 BC, 
30 BCo 
31 BC3 
32 BCy 
55 BCs 
34 BCg 
35 5 BMD ; 
36 6 BMD 9 
37 if BMD 3 
38 8 TC, 
39 TCo 
40 TC3 
41 Ty, 
42 TGs 
43 TC. 
44 9 TMD; 
45 10 TMD» 
46 11 TMD 3 
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In(Total Ingrowth) = bo + b,1n(BACL}) + bo1n(BACL>) cL b3BA!-° + 


byA,BA!+> + In(S) (6) 
where 
Bo = =7. 15817566 
b, = 0.706329932 
by = 1.97156496 


b3 = -8.04453668E-03 

by = -2.69619979E-03 

BACL, = basal area in the 4- through 6-inch diameter classes 
BACL» = basal area in the 7+-inch diameter classes 


BA total) basal area 


Ay time-since-last-cutting transform 


S = Minor's site index. 


The skewness statistic for this model was 1.10468, and the kurtosis statistic was 3.55524. 
These statistics indicate a moderate amount of "nonnormality.'' Examination shows the 
residuals were divided into two groups. One group represented the zero valued ingrowth 
observations (which has been set to a very small positive number when fitting the log models). 


The strange behavior of the residuals did not appear reasonable. Therefore, the anti- 
log of model 6 was fitted to the total ingrowth data through linear, least squares regression 
through the origin, to produce the following model: 


afPo + b3BA!-° + b,A,BA!-°) 


Total ingrowth = a9S(BACL,)”!(BACL,)?2 (7) 


where 

ag = 1.41650289 

bo, Dy, =--, Dy = Yegression coefficients from model (6). 
The skewness statistic for this model was 0.818989, and the kurtosis statistic was 2.44519. 
These values indicate that the residuals about model (7) are more normally distributed than 
those about model (6). Also, the residuals were in one group about the model, and they 
increased, indicating a necessity to weight the model. 

To determine the proper weighting scheme, the residuals were divided into seven predicted 


total ingrowth classes, and the variance for each class was then computed. Using these data, 
the following model was then developed: 


0.804660585 : 
(8) 


Variance = 8.42590476 + 3.11147074 (YHAT) 
where 


YHAT = predicted ingrowth from model (7). 


The reciprocal of predicted variance then provides the weights necessary in weighted regression. 
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Because of the behavior of the residuals, it was concluded that a weighted nonlinear 
regression model was the most appropriate model for total ingrowth. Using the model form 
and initial parameter values of model (7) and the weights derived from model (8), a weighted, 
nonlinear regression run was made to provide the following model: 


a(eaBA + csA,BA‘*) 


Total ingrowth = cgS{BACL,)°! (BACL2)°2 (9) 
where 

Co = 0.42994711 

C1, = 0.66118156 

Co = -0.46026739 

cz = -4.8091408E-06 

cy = 2.7316853 

Cs = -1.6118274E-06. 


This model behaves well except when BACL) approaches zero, which causes the model to 
"blow up.'' To control this, a value of 1.0 was added to BACL,, and the final weighted, non- 
linear model was fitted. 


Through-Growth Model Development 


Program RISK was used to develop the through-growth model. It was hypothesized that 
through-growth would most likely be correlated with predicted basal area growth in the 4-inch 
diameter class (BJBAGI1-4" and BJBAG2-4'') and with predicted total ingrowth (TINGRO). Using 
these two independent variables and the time-since-cutting variable, A,, RISK runs 1 through 
10 listed in table 36 were made. Analysis of these runs indicated that the two best independ- 
ent variables for predicting through-growth were BJBAG1-4"' (or BJBAG2-4") and A, *TINGRO. 
Using these independent variables, runs 11 and 12 were made, and the chi-square ''goodness-of- 
fit" statistics across predicted basal area growth classes and across predicted through-growth 
classes were computed. Across predicted basal area growth classes, the tests were insignif- 
icant at the 99 percent testing level (tables 37 and 38). Across predicted through-growth 
classes, the model with BJBAG1-4" was also insignificant; however, the model with BJBAG2-4" 
was Significant. This latter test result was very surprising because, up to this point, the 
two blackjack pine basal area growth equations behaved consistently in relation to each other. 
Comparison of the mean squared errors, "'t"'-values, regression coefficients, and mean predicted 
basal area growths all indicate that the two independent variables, BJBAG1-4'" and BJBAG2-4", 
and the two models do not differ greatly. This conclusion is further substantiated by the 
Similarity of the results for chi-square tests across predicted basal area growth classes. 

It was decided, therefore, to ignore this chi-square test and accept model 12 as being the 
final model for BJBAG2-4". 
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Table 36.--Listing of through-growth models tested 


eee Independent variables 

i! BJBAG1-4'"', TINGRO, A,*BJBAG1-4", A,*TINGRO, Aj, 
2 BJBAG2-4", TINGRO, A,*BJBAG2-4", A,*TINGRO, Aj 
3 BJBAG1-4'', TINGRO, Aj 

4 BJBAG2-4'', TINGRO, Aj 

5 BJBAG1-4'', TINGRO, A, *BJBAG1-4" 

6 BJBAG2-4'', TINGRO, A, *BJBAG2-4" 

7 BJBAG1-4'', TINGRO, A,*TINGRO 

8 BJBAG2-4"', TINGRO, A, *TINGRO 

9 BJBAG1-4'', TINGRO 
10 BJBAG2-4'"', TINGRO 
11 BJBAG1-4'', A,*TINGRO 
12 BJBAG2-4'', A, *TINGRO 


Table 38.--Chi-square test across predicted basal area growth classes for final through-growth 
model involving BJBAG1-4" 


Hee Total number Actual Predicted Chi-square 

asal area of in growth through- through- Baa 
growth class trees in class growth growth 

0.000 - 0.035 80 18 92 Sail 
2035; = .040 555 82 94.57 ib 
-040 - .045 715 134 139515 7) 
.045 - .050 493 12 105.78 4 
-050 = .055 614 133 141.42 =D 
-055 - .060 431 120 110.95 Sy 
.060 - .065 267 67 T0937 1 
.065 - .070 142 34 44.75 2.6 
.070+ 50 20 ZL 50S oak 


Chi-square statistic = 10.4 
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Table 38.--Chi-square test across predicted basal area growth classes for final through-growth 
model involving BJBAG2-4" 


ee Total number Actual Predicted Garson 

asal area of in growth through- through- 

growth class trees in class growth growth value 

0.000 - 0.035 80 18 AS At Bes) 
-03557= 040 555 82 95.19 ac 
.040 - .045 IAS 134 139.54 a2 
1045 °= 5050 493 ri2 105.75 .4 
2050 = 3055 614 133 ANS 2 55) 
A055: = .»..060 431 120 110.62 .8 
Z060) =") X065 267 67 75.63 a) 
065 =» 1.070 142 34 44.61 Zi 
070+ 50 20 21.86 a2 


Chi-square statistic = 10.3 
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APPENDIX H 


DEVELOPMENT OF HEIGHT EQUATIONS 
Even-aged Height Equation 


An examination of the Taylor Woods height data showed that average diameter class height 
increased asymptotically to the height predicted by Minor's (1964) height equation, as the 
diameter class size approached the maximum for the stand. Functionally, this behavior can be 
expressed: 


Dan 
Hp = MH - bj (MH-4.5)(1.0 - De 


With this function, average height of the diameter class (H,) approaches the height predicted 
by Minor's (1964) equation (MH) as diameter class size (D) approaches the maximum diameter 
class size that the stand has achieved in its development (DM). Conversely, H, approaches 
4.5*b, as D approaches zero. How quickly H approaches H, depends upon the value of n; the 

: : . Ae aD) 
larger its value the quicker H will approach Hp as D increases. 


The following model was used to determine the appropriate values of b, and n: 


H,-4.5 n 


D a D 
= b*(1.0 - =) 


eos (aoa 5 


By varying n between 1.0 and 2.5 at increments of 0.05, 30 different independent variables were 
formed. An all-combination screening run using least squares regression through the origin was 
then used to first select the best value of n and then to compute the regression parameter(b)). 
The resulting RMSQR for the final model was 0.1419. 


Uneven-Aged Height Equation 


A number of model forms have been proposed for characterizing height as a function of 
diameter, including: 


Seale Ome 


‘ b3(1nD)? 
4.5 + by*D 2+e (2) 


H 


4.5 + aj°e 


H 


The parameters K and m of model (1) have either been fixed at values of zero and one respec- 
tively (Curtis 1967; Embry and Gottfried 1971; Burkhart and Strub 1974), or they have been 
included as regression parameters (Monserud 1975). Values of the parameters have been deter- 
mined either through the logarithm transformation and linear regression process for models (1) 
and (2) (Curtis 1967; Embry and Gottfried 1971; Burkhart and Strub 1974), or they have been 
determined through nonlinear regression for model (1) (Monserud 1975). 


To incorporate the expected effect of site index upon these height-diameter models, the 
following models were designed: 


c3(D+K) ™ 
H = 4.5 + c,°S°2-e (3) 


Bye . Fie 
H = 4.5 + d,-$42-e9°*S* (D+K) (4) 


D, 
45 2 3-52 see 
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Models (3) and (4) are variations of model (1). Model (3) scales the basic height-diameter 
curves by the stand's site index, while model (4) also allows a site index related shift in 
the form of the height-diameter curve. Model (5) is a form of model (2) in which the basic 
height-diameter curve is scaled by site index. In applying these models to diameter classes, 
the dependent variable was defined as the average height of the Dth diameter class for a stand 
with site index of S. 


To fit models (3), (4), and (5) to the uneven-aged height data, they were first transformed 
by use of natural logarithms to the following: 


In(Hy-4.5) = f, + f2*InS + £3+(D+K) Ba, (6) 

In(H)-4.5) = 91 + g2inS + g3*S+(D+K) (7) 
2 

In(H)-4.5) = hy, + holnS + h3inD + hi(1nS) (8) 


For models (6) and (7), screening runs were then made with values of K ranging from 0.0 to 
40.0 and values of m ranging from 0.1 to 4.0. Model (8) was fitted using a single least 
Squares regression run. The three models were then tabulated over D and S to check for reason- 
ableness of behavior, their residuals were plotted over predicted average height to determine 
homogeneity of variance, and the residuals were also examined for normality. The final model 
form selected was (6) with K = 35.0 and n = 2.0. The residuals about this model were normally 
distributed with homogeneous variance, and therefore the lognormal bias correction proposed 

by Oldham (1965) and Baskerville (1972) could be used. This correction consisted of adding 
one-half the mean-squared error (0.030203189) to the intercept of model (6) and then taking 
the antilog of model (6) to obtain the final uneven-aged height model. As an index of fit, 
the RMSQR for the final log model was 0.1270. 
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APPENDIX | 
EXAMINATION OF STATISTICAL TESTS USED IN VALIDATION 


To pursue possible statistical procedures for testing differences between actual and 
predicted diameter distributions, I conducted a literature review of the validation techniques 
used in operations research. Of those methods described by Mihram (1972a, 1972b), Van Horn 
(1971), and Kleijnen (1974), two seemed appropriate for additional evaluation--the chi-square 
"goodness-of-fit" test and ordinary least squares regression. 


The chi-square "goodness-of-fit" test statistic can be computed as: 


where 
x2 = computed test statistic 
n = total number of observations 
Ny = number of observations that fell into the ith class 
p. = theoretical probablity of an observation falling into the ith class 


k = total number of classes. 


This statistic is normally used to test the frequency distribution of test data against a 
theoretical distribution (Snedecor and Cochran 1967; Kmenta 1971). Therefore, it seemed to be 
a natural statistic for testing an observed diameter distribution against a theoretical (or 
predicted) diameter distribution. An examination of the assumptions underlying this test, 
however, revealed problems that render the usage of the chi-square "goodness-of-fit" statistic 
for testing diameter distributions questionable, if not inappropriate. 


The foundation of the chi-square test is the multinomial distribution. If the assumptions 
of multinomial distribution can be met, then the limiting distribution of the chi-square 
"goodness-of-fit" statistic is the chi-square distribution. The two assumptions of the multi- 
nomial distribution of interest are: (1) the probability of an observation being in the ith 
class, p., is fixed; and (2) the number of observed values in each class must sum to n (Mood 
and others 1974; Kendall and Stuart 1973). 


Translating these assumptions to the problem at hand, the first assumption would state 
that the probability of a tree being in a specified diameter class is fixed. While this 
assumption appears to be unreasonable, it may be possible to define the p,'s as being condi- 


tional upon site and stand conditions. Therefore, p. would be the probability of a tree being 
in the ith diameter class at the end of the growth period given the present diameter distri- 
bution and the specified site index and time-since-last-cutting. 


The second assumption states that the total number of trees in the observed diameter 
distribution must equal the total number of trees in the predicted diameter distribution. To 
achieve this condition, upgrowth, ingrowth, mortality, and conversion from blackjack to yellow 
pine must be predicted exactly--and this is highly unlikely. Therefore, while the problem 
with the first assumption may be avoidable, violation of the second assumption seems unavoid- 
able. As a result, the use of the chi-square "goodness-of-fit" statistic for testing differences 
between predicted and actual diameter distributions is judged to be inappropriate. 
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The basic model for the ordinary least squares regression approach is: 


NGS Weer ap oN Oo sae Sc 
1 AL, Jt 
where 
Y, = actual value of a variable (or actual output of the model) 
Y. = predicted value of a variable (or predicted output of the model) 


a,b = regression coefficients 
e, = a random deviate. 


An intercept (a) of zero and a slope (b) of one would indicate, presumably, that the predicted 
output of the model is an unbiased estimate of the actual output of the model (Cohen and Cyert 
1961; Kleijnen 1974). If the random deviates are distributed normally with mean zero and 
hemogeneous variance, then standard tests can be used to determine if the intercept is signifi- 
cantly different from zero, and if the slope is significantly different from one (Draper and 


Smith 1966). Aigner (1972) has shown that the predictive model Y must be deterministic for 
this approach to be legitimate; and the simulator developed in this study is deterministic. 


In applying this technique to actual and predicted number of trees in a diameter class, 
the t-tests for determining the significance of the intercept and slope from zero and one, 
respectively, did not prove to be very useful for making comparisons and decisions. Tests 
were found in which the intercept and slope were greatly different from their expected values, 
but, because of the high degree of impreciseness between actual and predicted values (that is, 
because of the high MSE of the regression), the tests were "'insignificant.' Conversely, an 
intercept and/or slope not greatly different from its expected value was often found to be 
"significantly" different if a low MSE existed (that is, where predictions were highly precise). 
This limitation of the regression t-tests is the same as attributed by Freese (1960) to the 
paired t-test. 


Another statistical test considered for use in validation was the Kolmogorov-Smirnov (K- 
S) test. The standard K-S test was developed for continuous distributions. Because the 
diameter distribution data are discrete, application of the standard K-S test to it would have 
resulted in conservative tests (Conover 1971). A discrete K-S test had been reported in the 
literature (Conover 1972; Horn 1977), but the time and expense of developing and applying the 
required computer program was judged too great for this study. 
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APPENDIX J 
DESCRIPTION OF CONTROL CARDS 


Field function Value Column Name Format 
Card 1 

Title of job 1-80 TITLE 20A4 
Card 2 

Minor's (1964) site index XXX. X 1-5 S BSe tt 


Number of 5-year growth periods 
in run (<80) KX 7-8 NGP 12 


Number of 5-year growth periods since 
last cutting in stands (1<TCUT<60). 


Set to 60 for virgin stands. XX 10-11 TCUT 12 
Breast height stand age. Set to zero 
for uneven-aged stands. XKX 13-15 AGE F3.0 
Cutting scheme 
No cutting 0 17 JCUT Il 
Removals are standard for all cuts and 
are designated by diameter classes. i 
Removals are specified for each cut and 
are designated by diameter classes. 2 
Removals are specified for each cut and 
are designated by thinning type. 3 


Cutting specifications 
Residual trees per diameter class is 
specified in number of trees 
(to be used with JCUT=1 or 2). 1 19 SCUT I] 
Residual trees per diameter class 
is specified as a proportion 


(to be used with JCUT=1 or 2). 2 
Residual total stand basal area is 

specified (to be used with JCUT=3). 55 
Residual total number of trees in stand 

is specified (to be used with JCUT=3). 4 


Proportion of mortality to be recovered 
in thinnings VERIO O OS 21-26 MORREC F6.4 


Proportion of removals to come from below 
when thinning from both above and below 
(to be used when JCUT=3 and thinning 
type =4). Set to 0.0 if JCUT#3 or 
thinning type #4. X.XXXX 28-33 PROREM F6.4 
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A "C'' punched in columns whose ordinals 
are the growth periods for which cutting 
to occur. Blank if no cutting is to 


is 


Field function 


Cards 


occur. 


CGandiGs)M4s eae SNL 


From 1 to as many as 6 cards 


Initial blackjack pine array 


Note: Omit if no blackjack pine. 


For each card: 


Card 


Card type 

Diameter class (<45) 

Number of trees in diameter class 
Diameter class (<45) 

Number of trees in diameter class 


Diameter class (<45) 
Number of trees in diameter class 


Card Nl+l 


type 
(This card must be present, 
even if no blackjack pine) 


Gard (s)i NUP2yi1.2 NZ 


From 1 to as many as 6 cards 


Initial yellow pine array 


Note: Omit if no yellow pine 


For each card: 


Card 


Card type 
(Rest of card same as for 


initial blackjack pine cards) 


Card N2+1 


type 
(This card must be present, 
even if no yellow pine) 


Gandi(s)i N2Z225 oe NS 


Value 


C or blank 


INITIAL BJ 
XX 

XXX. XX 

XX 

XXX. XX 


XXX. XX 


INITIAL YP 


DONE 


For JCUT=1 or 2, only (for JCUT=3, go to 


next definition of these cards) 
From 1 to as many as 6 cards 


Blackjack pine cutting array 


Note: Omit if no blackjack pine cutting 
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Column 


1-80 


1-10 
12-13 
15-20 
12-23 
25-30 


72-73 
75-80 


1-4 


1-10 


1-4 


Name 


NC 


CHECK 
ID(1) 
GD) 
ID(2) 
23) 


1D(7) 
T(7) 


CHECK 


CHECK 


CHECK 


Format 


80A1 


A4, 6X 
12 
ROg 2 
12 
EGR 


12 
F6.2 


A4 


A4, 6X 


A4 


Field function Value Column Name Format 


For each card: 


Card type CUTTING BJ 1-10 CHECK A4, 6X 
Diameter class (<45) XX 12-13 ID(1) 12 
If SCUT=1, number of residual 

trees in diameter class XOXO PAEX: 15-20 T(1) F6.2 
If SCUT=2, proportion of diameter 

class to be uncut X. XXXX 15-20 Gls) F6.4 
Diameter class (<45) XX 22-23 ID(2) 12 
If SCUT=1, number of residual 

trees in diameter class XXX. XX 25-30 1) EG@).2 
If SCUT=2, proportion of diameter 

class to be uncut SORKIN 25-30 T(2) F6.4 
Diameter class (<45) Ge 72-73 ID(7) I2 


If SCUT=1, number of residual 
trees in diameter class XXX: 75-80 Th (Gi) F6.2 


If SCUT=2, proportion of diameter 
class to be uncut » 9, 0,0,0.4 75-80 TGA) F6.4 


Repeat the above card type until 
all diameter classes are defined. 


Cardi(S)PNZ52e DNS 
For JCul—se only (for JCUT 1 or 2, go 
to previous definition of these cards) 
From 1 to as many as 27 cards 


Thinning array 
Note: Omit if no thinning 


For each card: 
Card type THINNINGS 1-9 CHECK A4, 5X 
Growth period for thinning XX 12-13 ID(1) 12 
(Order of thinnings on card 
need not be sequential) 


Thinning type: 


Thinning from below’ 1 15 ID(2) Tl 
Thinning from above® 2 
Thinning from below and above 3 
Thinning proportionally across 
all diameter classes? 4 


?Thinning from below is the removal of trees from the smallest diameter class first. 

8Thinning from above is the removal of trees from the largest diameter class first. 

9Proportional or random thinning is the removal of basal area (SUC=3) or trees (SCUT=4) 
proportionally from all diameter classes. The form of the basal area or diameter distribution 
after thinning will be the same as that before thinning. 


Field function 


If SCUT=3, residual blackjack pine 
Basal area after thinning 


If SCUT=4, residual blackjack pine 
Number of trees after thinning 


If SCUT=3, residual yellow pine 
Basal area after thinning 


If SCUT=4, residual yellow pine 
Number of trees after thinning 


Growth period for thinning 
Thinning type 


If SCUT=3, residual blackjack pine 
Basal area after thinning 


If SCUT=4, residual blackjack pine 
Number of trees after thinning 


If SCUT=3, residual yellow pine 
Basal area after thinning 


If SCUT=4, residual yellow pine 
Number of trees after thinning 


Growth period for thinning 
Thinning type 


If SCUT=3, residual blackjack pine 
Basal area after thinning 


If SCUT=4, residual blackjack pine 
Number of trees after thinning 


If SCUT=3, residual yellow pine 
Basal area after thinning 


If SCUT=4, residual yellow pine 
Number of trees after thinning 


Repeat the above card type until all 


thinnings specified on card 3 are defined. 


Value 


XXX .XX 


XXXX.X 


XXX. XX 


XXXX.X 


XXX. XX 


XXXX .X 


XXXX.X 


XXX. XX 


XXXX.X 


XXX. XX 


XXXX.X 
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Column 


18-23 


18-23 


38-43 


45-50 


45-50 


52-53 


55 


58-63 


58-63 


65-70 


65-70 


Name 


T(1) 


TL) 


T(2) 


11) 
ID(3) 


ID(4) 


T(3) 


T(3) 


T(4) 


T(4) 
ID(5) 


1D(6) 


T(S) 


T(S) 


T(6) 


T(6) 


Format 


F6.2 


Fo.1 


Field function Value 
Card N3+1 


Card type DONE 
(This card must be present, 
even if no blackjack cutting or thinning) 


Gard(s) N3+2;,....,N4 
For JCUT=1 or 2 only 
From 1 to as many as 6 cards 


Yellow pine cutting array 
Note: Omit if JCUT=3 or if no 
yellow pine cutting. 


For each card: 
Card type CUTTING YP 
(Rest of card same as for 


blackjack pine cutting cards) 


Repeat the above card type until 
all diameter classes are defined. 


Card N4+1 


Card type DONE 


(This card must be present 
if JCUT=1 or 2) 


If JCUT=2, then cards (N2+2) to (N4+1) 
are repeated for as many times as 
ehere are “G"s" punched in card 3. 


Card N4+2 


Card type if another growth 
analysis is to be done CONT INUE 


Card type if last growth analysis FINISH 


Column 


1-4 


1-10 


1-4 


1-8 


1-6 


If CONTINUE, then repeat cards 1 to (N4+1) for the new analysis. 


Ww 
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Name 


CHECK 


CHECK 


CHECK 


CHECK 


CHECK 


Format 


A4 


A4, 6X 


A4 


A4 


A4 
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